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CHAPTER  1 
INTRODUCTION 

With  the  increasing  use  of  electronic  computers,  random  sampling 
methods  have  become  a  very  useful  tool  for  providing  solutions  to  prob¬ 
lems  involving  probability  as  well  as  for  problems  of  a  deterministic 
nature.  The  availability  of  sequences  of  numbers  which  appear  to  be 
drawn  at  random  from  particular  probability  distributions  is  a  vital 
ingredient  in  the  random  sampling  process.  Such  numbers  are  referred 
to  as  pseudo-random  numbers  or,  for  convenience,  simply  as  random  num¬ 
bers.  This  paper  is  concerned  with  the  rapid  and  accurate  generation 
of  such  numbers  in  a  stored  computer  program. 

The  use  of  random  sampling  to  estimate  distribution  functions  origi¬ 
nated  with  Student  [17]  in  1908.  At  that  time  and  for  some  time  to 
follow,  necessary  random  numbers  were  obtained  by  drawing  cards  from  a  deck, 
counters  from  an  urn,  or  by  rolling  dice.  Such  processes  are  very 
slow  and  make  it  quite  difficult  to  insure  randomness. 

To  facilitate  the  use  of  random  numbers  large  tables  of  random  digits 
were  compiled.  The  first  such  table,  published  by  Tippett  [19]  in  1927, 
consisted  of  41,600  digits  taken  at  random  from  census  reports  and  com¬ 
bined  into  10,400  four  digit  numbers.  The  requirement  for  a  larger  set 
of  numbers  led  to  the  publication  of  100,000  random  digits  in  1939  by 
Kendall  and  Babington-Smith  [7].  These  digits  were  generated  by  means  of 
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a  mechanical  process  and  were  the  first  to  be  so  produced.  Kendall  and 
Babington-Smith  are  also  responsible  for  developing  many  of  the  tests 
frequently  applied  to  sequences  of  random  numbers  (5,6).  Other  such 
tables  have  since  been  constructed,  primarily  by  the  use  of  physical 
devices,  culminating  in  the  most  extensive,  "A  Million  Random  Digits 
with  100,000  Normal  Deviates"  published  ty.the  RAND  Corporation  [15]  in 
1955. 

The  use  of  a  physical  device  for  the  generation  of  random  numbers 
on  line  with  a  computer  is  both  expensive  and  difficult  to  maintain.  The 
storing  of  a  table  of  random  numbers  on  magnetic  tape  or  on  cards  is  also 
an  unsatisfactory  method  of  generating  numbers  for  use  in  a  computer. 

This  would  necessitate  the  use  of  an  input  device;  also  the  time  required 
to  read  in  the  numbers  would  be  excessive. 

The  development  of  arithmetic  procedures  for  the  generation  of  ran¬ 
dom  numbers  began  in  the  1940's  with  the  introduction  of  computers.  The 
first  such  procedure  was  suggested  by  von  Neumann  and  Metropolis  in  1946 
and  is  described  in  [4],  It  was  also  during  this  time  that  calculations 
involving  random  numbers  received  the  picturesque  name  "Monte  Carlo."  The 
"middle-square"  method  proposed  by  von  Neumann  and  Metropolis  is  simple, 
fast,  and  requires  only  an  initial  starting  value.  In  this  procedure  each 
new  number  is  produced  by  taking  the  middle  n  digits  of  the  square  of  the 
previous  n-digit  number.  As  pointed  out  in  [4]  and  in  [16]  the  sequence 
of  numbers  produced  using  this  method  sooner  or  later  degenerates  to  a  cy¬ 
cle  which  often  is  very  small,  and  at  worst  consists  of  only  a  single  num¬ 
ber.  In  addition  some  of  the  statistical  tests  performed  on  samples  gen¬ 
erated  by  the  middle-square  procedure  have  resulted  in  failures. 
Improvements  on  the  von  Neumann  method  are  plagued  by  similar  difficulties 
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[2]  and  [20]. 

Of  the  various  arithmetic  procedures  sequences  of  random  numbers 
with  the  best  statistical  properties  and  longest  periods  are  generated 
.  by  means  of  the  congruence  relation 

=  a  •  +  c  *(mod  m) .  (1.1) 

The  representation  (1.1)  is  termed  the  mixed  congruential  method.  The 
multiplicative  congruential  method  is  defined  by  taking  c  =  0  in  (1.1) 
Lehmer  [8]  is  credited  with  the  invention  of  the  multiplicative  method. 

In  their  very  informative  survey  Hull  and  Dobell  [3]  prescribe  conditions 
for  a,  m,  c,  and  xq  which  insure  maximum  period.  Other  procedures  for 
the  genration  of  random  numbers  based  upon  reduced  Fibonacci  series  [21] 
and  upon  transcendental  numbers  [3]  are  inferior  to  those  based  on  the 
congruence  relation  (1.1)  (See  [4]).  In  [9]  MacLaren  and  Marsaglia 
present  a  table  look-up  scheme  which  seems  tQ  offer  some  promise. 

Random  numbers  generated  by  arithmetic  procedures  are  not  "truly 
random"  in  that  the  entire  sequence  is  determined  in  advance  and  can 
be  reproduced  simply  by  using  the  same  starting  value,  x0.  The 
statistical  behavior  of  sequences  generated  by  both  the  mixed  and 
foultiplicative  congruential  methods  is  quite  good,  and,  in  fact, 
numbers  genrated  in  this  fashion  do  appear  to  be  drawn  at  random  from 
the  unifoi'm  (0,1)  distribution 

f(u)  =  1,  0<u<l .  (1.2) 

The  reproducibility  of  scqucnc.es  is  an  advantage  in  debugging  and  in 
certain  calculations  it  may  be  desirable  to  reproduce  a  given  sequence. 
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Bargmann  [1]  has  provided  a  procedure  for  the  generation  of  indepen¬ 
dent  uniform  (0,1)  random  numbers  on  a  binary  computer.  This  proce¬ 
dure,  which  is  a  form  of  the  multiplicative  congruential  method,  is  defin¬ 
ed  by  .  ' 

u  =  a  •  u  (mod  232)  (1.3) 

n+1  n v  J  v  J 

where  a  is  chosen  such  that  a  +  1  and  a  -  1  end  in  as  few  zero  bits  as 

possible.  This  is  insured  by  requiring  that  neither  a  -  1  nor  a  +  1  be 

divisible  by  2  for  k  =  3,4,....  Choosing  a  as  an  odd  power  of  5 

determines  that  the  low  order  3  digits  will  be  125,  and  neither  124  nor 

k  . _ 

126  is  divisible  by  2  for  k  =  3,4, _  Hence,  a  =  5 1 i  is  a  reasonable 

choice.  This  procedure  r*q*"^es  that  uq,  the  starting  value,  be  an 
odd  positive  integer. 

The  use  of  this  procedure  results  in  a  very  fast  computer  program 
requiring  only  5  operations:  load,  multiply,  and  store.  A  description 
of  the  usage  of  this  program  along  with  time  and  storage  requirements 
is  given  in  Chapter  6  of  this  paper.  Results  of  statistical  tests 
performed  on  this  procedure  are  presented  in  Chapter  5,  and  a  listing 
of  the  program  is  given  in  Appendix  A. 

The  preceding  discussion  has  dealt  only  with  the. generation  of 
random  numbers  which  appear  to  be  uniformly  distributed.  In  principle 
it  should  be  very  easy  to  obtain  any  other  distribution  from  the  uni¬ 
form  distribution,  requiring  only  a  solution  to  the  equation 

u  =  F(y)  (1.4) 


for  y,  where  u  is  uniformly  distributed  on  the  interval  (0,1)  and  F  is  the 
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required  cumulative  distribution  function.  When  the  inverse  of  F  is 
known,  as  for  the  exponential  distribution,  this  is  a  simple  matter: 

y  =-ln(u) .  (1.5) 

However,  the  evaluation  of  the  In  function  may  be  somewhat  time  consuming. 
An  alternative  approach  might  be  to  store  a  value  of  y  for  each  possible 
value  of  u  based  on  the  relation 

y  =  F_1(u).  (1.6) 

To  generate  y,  simply  generate  a  uniform  (0,1)  random  number  u  which  will 
determine  the  location  of  a  stored  value  of  y.  Let  u  be  given  to  8 
digits  and  let  8  be  the  base  of  the  number  system  in  which  the  digits 
are  represented,  then  for  this  procedure  a  total  of  08  storage  locations 
are  required. 

Various  approximations  have  been  proposed  for  the  normal  distribu¬ 
tion.  Most  of  these  involve  taking  the  sum  of  a  fixed  number  of  uniform 
(0,1)  deviates  or  tie  use  of  Chebyshev  approximations  [14]  or  a  combi¬ 
nation  of  the  two  [18].  Such  approximations  frequently  lack  accuracy 
and  are  either  slow  or  space  consuming. 

The  procedures  suggested  by  Marsaglia  ct  al.  in  [11], [12],  and  [13] 
for  transforming  uniform  (0,1)  random  numbers  to  random  numbers 
having  other  distributions  are  superior  to  3ny  encountered.  They  are 
fast,  require  minimal  space,  and  arc  simple  to  program.  In  addition 
these  procedures  are  completely  accurate',  the  precision  of  the  result 
is  dependent  only  on  the  word  size  of  the  computer.  These  procedures 
along  with  programs  and  results  of  statistical  tests  arc  presented  in 
detail  in  this  paper. 


CHAPTER  2 


'GENERATION  OF  DISCRETE  RANDOM  VARIABLES 
IN'  A  COMPUTER 

Let  Y  be  a  discrete  random  variable  with  point  probabilities 
=  Pr(Y  ~  y.)  for  i  -  1,2,....  The  direct  way  to  generate  Y  in  a 
computer  is  to  generate  a  uniform  (0,1)  random  number  u  and  put  Y  -  y^ 
if  pi  +  p?  +  ...  +  p.  ^  <  u  <_  pi  +  p2  +  ...  +  .  However,  techniques 

based  on  this  method  lead  to  complicated  programs  that  are  excessively 
time  consuming. 

An  alternative  method  proposed  by  Marsaglia  [11]  is  simple  to 
program  and  requires  minimal  time  and  storage.  Let  for  i  =  1,2,..., 
n  be  expressed  by  k  digits  as  p.  =  - <5 1 ^ 62 i  •••  ^  the  domain  of 

the  random  variable  is  infinite,  the  probability  distribution  must  be 
truncated  at  some  p^.  The  fastest  method  for  generating  Y  is  as  follows. 

Let  6  be  the  base  of  the  number  system  in  which  the  <$^'s  are  repre- 

}- 

sented.  In  memory  locations  0  to  (3  -1)  store  ...  <5]<i  yj's, 

6j2<522  •  •  •  5k2  '  i ,  iin^n  •  •  •  5kn  Yn's-  If  u  is  a  uniform  (0,1) 

random  number,  u  -  . d jd 2  dg,  look  up  the  number  in  location  djd2...  dj. 
and  lot  that  be  Y.  j 

Though  this  may  be  the  fastest  method,  it  clearly  requires  an  exes- 
sive  amount  of  storage  space.  Even  if  the  pVs  are  truncated  to  four 
digits,  S1'  memory  locations  will  be  required. 
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Marsaglia  [11]  suggests  a  technique  that  offers  a  considerable  re¬ 
duction  in  storage  space  with  very  little  sacrifice  in  execution  time. 
For  convenience  assume  that  the  p^'s  are  truncated  to  four  digits.  We 


have  the  following  situation: 

Value  of  V 


Frobabi lity 


6  6  6  6 
3  J  2  )  3141 


6  6  6  6 
12  22  32  42 


6,  6„  6  . 

In  2n  3n  4n 


Define 


Fq  =  0, 


-r 


n 

£  6  . 

.  ,  n 

i=l 


(2.0.1) 


for  r  =  1 ,2,3,4,  and 


s 

Z 


n 

Z  6.. 


n  = 
s  j=i  i=i 


(2.0.2) 


for  s  =  1, 2,3,4  and  n0  =  0.  Segment  memory  locations  0  to  (JI4-1) 

into  four  mutually  exclusive  sets  such  that  set  1  consists  of  locations 

0  to  (Hj-1),  set  2  comprises  locations  Jl i  to  (Jl^-l),  and  in  general 

the  st^1  set  occupies  memory  locations  Rs  ^  to  (n^ - 1 )  for  s  =  1,2, 3,4. 

In  each  set  s  for  s  =  1,2, 3, 4  each  y.  is  stored  in  5  .  locations  for 

1  si 

i  =  1,2,  ...,  n,  a  total  of  II s  -  ns_^  locations.  The  total  memory  rc- 

4  n 

quirement  is  then  n4  =  E  E  5..  locations.  Mow  choose  set  1  with 

j=l  i=l 

probability  Pj,  set  2  with  probability  Pn, 


set  3  with  probability  P  -L ,  or 


set  4  with  probability  P,4.  Having  chosen  a  set,  select  at  random  a  loca¬ 
tion  within  that  set;  the  number  occupying  that  location  is  the  desired 
random  variable. 

This  procedure  gives  the  required  discrete  distribution,  as  can  be 
seen  by  defining; 

* 

n 


for  i  -  1,2,  ....  n  and  r  =  1, 2,3,4.  Thus  is  simply  the  probability 
of  choosing  y^  from  set  r.  Then  the  probability  of  generating  V  =  y.  is 


4 

£ 

r=l 


P 


r 


‘n 


4 

z  a"r<s  . 

n 


.6.6.66  . 

U  ?i  31  4l 


> 


which  is  the  probability  P(Y  =  yp. 

In  order  to  select  the  proper  set  and  location  within  that  set, 
generate  a  uniform  (0,1)  random  variable  u  =  .dp,,...,  and  let  A{  j } 
denote  the  contents  of  memory  location  j.  Then  if 

s-1  s 

E  P  <  u  <  E  P 
r=0  r  r=0  r 

put 

c  S_1 

Y  =  A  (djd  •  ••  d  +  n  -  0S  Z  p  ). 

r=0  T 


The  following  example  illustrates  this  procedure.  The  hexadecimal 
(6=16)  number  system  is  used  as  it  is  the  representation  of  the  IBM 
360/66, 


i 


1  0  O.AOOO 

2  '  1  0. 3800 

3  2  0.14A2 

4  3  0.OF5O 

5  4  0.0202 

6  s  0.020C 


Defining  and  Jig  as  before  gives 
P0  =  0 


Pi  =  16'1  E6  =  0.£ 
i=l11 

6 

P2  =  16‘2  16  .  =  (),  ip 
i=l21 

6 

P3  =  16"3  =  0, oof 

i=l  1 
’  6 

PH  =  lo'4  E6l)i  '=  0.0010 


n0  =  0 

1  6 

nl  =  r  F  6.  =  E 

.  j=l  i= 1 


and 


2  6 

jlo  =  E  E  6.  .  =  2D 


j=l  i=l 


3  6 

n,  =  I  I  6. ,  =  3C 


j=l  i=l 


4  6 

Jl4  =  E  E  s  4C. 

j=l  i=l 


J1 


The  four  sets  are  stored  as  follows: 


A  -  TABU; 
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SOT  1 

SET  2 

SET  3 

SET 

4 

Loc. 

Con. 

Loc. 

Con . 

Loc . 

Con. 

Loc. 

Con . 

Loc. 

Con. 

0 

0 

E 

i 

IE 

3 

2D 

2 

3C 

2 

1 

0 

F 

i 

IF 

3 

2E 

2 

3D 

2 

'  2 

0 

10 

i 

20 

3 

2F 

2 

3E 

4 

3 

0 

11 

i 

21 

3 

30 

2 

5F 

4 

4 

0 

12 

i 

22 

3 

31 

2 

40 

5 

5 

0 

13 

i 

23 

3 

32 

2 

41 

5 

6 

0 

14 

i 

24 

3 

33 

2 

rsi 

'cr 

5 

7 

0 

15 

i 

25 

3 

34 

2 

43 

5 

8 

0 

16 

2 

26 

3 

35 

2 

44 

5 

9 

0 

17 

2 

27  , 

3 

36 

2 

45 

S 

A 

1 

18 

2 

28 

3 

37 

3 

46 

5 

B 

1 

19 

2 

29 

4 

38 

3 

■47 

5 

C 

1 

1A 

3 

2A 

4 

39 

3 

48 

5 

D 

2 

IB 

3 

2B 

5 

3A 

3 

49 

5 

1C 

5 

2C 

5  - 

3B  , 

3 

4  A 

5 

.  - 

ID 

3 

.  .  ■  — 

4B 

5 

TABLE  2,1 
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To  generate  the  desired  random  variable,  let  u 
uniform  (0,1)  random  number  and  let  A{ j )  denote,  the 
location  j  and  proceed  as  follows: 


=  .did2  ...  be  a 
contents  of  memory 


1)  If  0  <  u  <  O.E  put  Y  =  A{ d 2 } 


2)  If  0.E0  <_  u  <  0.  FF  put  Y  *  A{d1d2  -  D?} 


3)  If  0.  FFO  <_  u  <  O.FFF  put  Y  -  A{d]d2d3  -  FC3) 

4)  If  O.FFFO  1  u  put  Y  =  A{d1d2d3d4  -  FFK4}. 


Examples : 


u  =  0.2170  ...  Y  =  A{ 2 }  =  0 

u  -  O.EFIO  ...  Y  *  A{EF-D2]  =  A{1D}  =  3 

u  =  0.FFFE  ...  Y  =  ACFFFE  -  FFB4}  =  a(4A)  =  -5. 


In  this  example  4C(]6)  memory  locations  are  required  as  compared  to 
16'*  memory  locations  for  the  fastest  method.  Suppose  the  times  for  cer¬ 
tain  operations  in  the  computer  are: 


Operation 


Time 


Compare  two  integers 

Subtract  two  integers 

Look  up  an  addressed  location 


Then  the  fastest  method  for  generating  V  requires  a  total  time  of  P  .  L 
and  10'.  memory  locations.  In  the  example  presented  only  4C(16)memory 
locations  are  required  and  an  average  generation  time  of 
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\ 


O.E  •  (P  +  L)  +  O.F  •  (2P  +  S  +  L)  +  0. OOF  •  (3P  +  S  +  L)  +  0.001 

•  (3P  +  S  +  L)  =  1.21,,..  •  P  +  0.20,  •  S  +  L 

llfjl  Uoj 

A  schemacic  for  storing  the  four  sets  of  discrete  variables  and 
the  subsequent  generation  of  the  desired  random  variables  is  presented  in 
figure  2. 1 . 
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FIGURII  2.1 

Flowchart  for  the  Generation  of  Random 
Numbers  from  a  Discrete  Distribution 

P  is  a  vector  of  probabilities  (p.  =  .6, .  . . .) 

Y  is  a  vector  of  the  corresponding  discrete*  variables  (integer  mode) 
n  is  the  number  of  elements  in  P  and  in  Y. 

$  is  the  base  of  the  number  system  in  which  the  6^ ^  ' s  are  represented 
(for  the  360/65  6  =  16,  for  a  decimal  machine  6  =  10). 

LGT  is  the  maximum  length  of  the  A-Table.  If  LGT  is  exceeded,  the  table 
loading  process  ceases  and  N,  is  adjusted  (since  N  was  calculated  on  the 
basis  that  loading  would  be  complete). 


GENERATE 


CHAPTER  3 


GENERAL  TIiCllNI  QUHS  1:0R  TUP.  GENERATION 
OP  PSEUDO- RANDOM  NUMBERS  HAVING  A 
CONTINUOUS  DISTRIBUTION  FUNCTION 

Certain  techniques  of  a  rather  general  nature  are  used  in' some  of 
the  routines  described  in  this  paper.  They  are  especially  well  suited 
for  use  in  very  fast  computer  programs.  The  principle  of  each  technique 
is  described  below  and  reference  will  be  made  to  these  techniques  as  they 
are  used  in  the  specific  generator  programs. 

3.1  The  Composition  Technique 

'  The  fundamental  procedure  used  for  generating  pseudo-random  numbers 
from  a  continuous  probability  distribution,  as  suggested  by  Marsaglia 
[10],  is  based  upon  a  decomposition  of  the  density  function  f  into  a 
mixture  of  3  densities: 

f(0  ~  piE;U)  +  r.;;:  (t)  +  P3E3OO  (3.1.1) 

where  pj  >  p2  ?  P3,  Pi  +  P2  +  P3  ~  1  :md  Pi  very  close  to  1. 

A  random  number  from  f(t)  is  obtained  as  either  a  number  from 
gl(t),  from  g2(0  or  from  gj(t)  with  respective  probabilities  pj,  p2,  p2. 
Consider  a  density  f(t)  defined  from  all  t  >_  a  (figure  3.1). 
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FIGURU  3.1 

The  density  gi(t)  is  represented  by  the  rectangles,  appropriately 
standardized,  including  the  shaded  upper  portions,  and  is  defined  for 
a  <  t  <  b.  The  width  of  the  rectangles  is  A,  a  quantity  whose  value 
is  dictated  by  the  number  system  used  within  the  computer.  For  a  binary 
machine  using  hexadecimal  arithmetic,  A  =  0.1^16^  =  0.0625 
The  total  area  represented  by  the  rectangles  is  pi  which  clearly  is 
very  close  to  1.  The  density  g 2  C t )  is  represented  by  the  "triangular" 
regions,  appropriately  standardized,  lying  below  f(t)  and  above  g)(t), 
and  is  also  defined  for  a  <_  t  <  b.  The  total  area  occupied  by  the 
"triangles"  is  p?.  The  density  g3(t)  is  defined  for  t  and  represents 
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the  tail  of  l(t).  Its  area  is  p3. 

Random  variables  having  density  gi(t)  may  be  rapidly  generated  an 
follows : 

Define 


and 


N  -  (b-a)/A 


=  a  +  (i-1)  •  A  for  i  =  1,2,3,  ...,  N. 


(3.1.2) 


(3.1.3) 


Now  assign  each  discrete  t^  a  probability 


P.  =  A  •  f (t.  ,)  =  .6,  .6  . . 
l  ^  i  +  lJ  l i  2i 


(3.1.4) 


The  Pi's  are  simply  the  areas  of  the  individual  rectangles  comprising 
gj(t).  Based  on  (say)  the  high  order  4  digits  of  the  P^'s  store  the 
discrete  t. 's  in  a  table  according  to  the  technique  described  in 
Chapter  2  of  this  paper.  Thus  to  generate  a.  random  variable  Y  from  this 
portion  of  gi  (unshaded  in  figure  3.1),  generate  a  uniform  (0,1)  random 
number  u  =  .d^  ....  If 

N 

u<  l  • 6  6  . 6, .6. 

i=1  U  21  31  41 

then  allow  the  4  high  order  digits  of  u,  did2d3d4,  to  locate  a  particular 
from  the  discrete  tabic,  as  described  in  Chapter  2,  and 
set  Y  =  ti  +  A  •  ( . d sd c, . . . ) . 

Random  variables  from  the  residual  of  g>,  corresponding  to  the  re¬ 
maining  digits  of  the  P^'s  (. 0000-5  ^6  ^ )  and  represented  by  the  shaded 
region  in  figure  3.1,  are  generated  according  to  the  following  scheme: 

Let  N,  t^,  and  be  defined  as  before  (3.1.2),  (3.1.3),  and  (3.1.41; 


and  define: 
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H.  -  .6  .6  .6  .6  . 
1  11  21  31  4l 


th—  — 

(area  of  the  unshaded  portion  of  the  i  rectangle), 


(3.1.5) 


R. 

l 


P. 

i 


H. 

l 


(area  of  the  shaded  portion  of  the 


.th 

l 


rectangle) , 


(3.1.6) 


Ti  *  FCti+.l)  •  F(V  '  pi  (3.1.7) 

where 

t. 

F(t . )  =  /  1  f(t  )dt. 

1  1 

(T\  is  the  area  of  the  "triangular"  region  above  the  i  rectangle),  and 

N 

P!  =  £  H.  (3.1.8) 

1  i=l  1 

(total  area  occupied  by  the  unshaded  rectangles).  It  is  obvious  from 
the  above  definitions  that 

N 

p2=  IT.  ,  (3.1.9) 

i=l  1 

N 

Pj  =  IP.  ,  and  (3.1, 10) 

i=l  1 


N 

pi  '  p!  =  .^Ri  (3.1.11) 


(Sce(3.1.1)  for  the  significance  of  pj  and  p2).  Now  store  the  t.^'s  in 
N  consecutive  locations.  The  order  in  which  the  t^'s  are  stored  is  not 
important  .  However,  for  convenience  let  D ( 1 )  denote  the  t-value  occupying 


Tlu-y  may  be  overlapped  with  the  previously  stored  t-Valucs  in  order 
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the  first  location,  D(2)  denote  the  t-value  occupying  the  second  loca¬ 
tion,  etc.  For  each  t^  define  a  pair  of  values  C(k)  and  B(k)  as  follows: 

C (1 )  =  pj  *  R[D(1)],  (3.1.12) 

B(k)  -  C(k)  +  T[D(k)J,  and  (3.1.13) 

C(k+1)  =  C(k)  +  RlD(k^l)]  (3.1.14) 

for  k  =  1,2,  N.  R [D (k) J  and  T ID Ck) 3  denote  respectively  the  area 

of  the  rectangle  residue  (shaded  in  figure  3.1)  and  the  area  of  the 
"triangular"  region  corresponding  to  the  t-value  occupying  location 
D(k).  It  is  apparent  that  . 

B(N)  =  Pj  +  p2.  (3.1.15) 

If  u  is  a  uniform  (0,1)  random  number  such  that  p|  f.u  <  pj  p2, 
a  random  variable  Y  is  generated  from  either  the  rectangular  residues 
of  gj  or  the  "triangular"  regions  of  g2  as  follows: 

If  B(k-l)  £u  <  C(k)  then  generate  a  new  uniform  (0,1)  random  number 
v  and  set  Y  =  D(k)  +  A  •  v.  A  random  variable  generated  in  this  fashion 
will  have  the  density  represented  by  the  shaded  rectangles  in  figure  3.1. 
If  C(k)  <_  u  <  B(k),  then  the  random  variable  Y  must  be  generated  from 
the  "triangular"  density  of  g2  corresponding  to  the  t-value  occupying 
location  D(k),  The  method  used  depends  upon  the  particular  function 
f(t).  The  same  applies  for  random  variables  from  g3  for  t  >  b. 
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The  techniques  used  for  the  exponential  and  normal  densities  will  be 
described  subsequently. 

Since  p]  is  close  to  1,  most  of  the  time  a  random  variable  from 
gi  is  generated.  Though  g 2  and  03  may  be  complicated  and  require  longer 
running  time,  the  average  time  per  generation  is  small  since  they  must 
be  handled  so  rarely. 

3.2  The  Acceptance-Rejection  Method 
Consider  the  density  h(t)  defined  on  the  interval  (bo»£i).  See 
figure  3.2,  Now  if  u,v  is  a  pair 

c 

h(t) 

y 


$0 

FIGURE  3.2 

of  independent  uniform  (0,1)  random  numbers,  x  =  £0  +  u  *  (£]-£0)  and 
y  “  v  •  c  will  have  uniform  densities  on  (£,o*Sl)  anc*  (0,c)  respectively. 
Suppose  that  y  <_h(x).  Then  the  conditional  distribution  function  of  x 
is  given  by 

x  h  (t )  x 

P[X  <  x  |  Y  <  hix)]  »  /  /  dydt  =  /  h(t)dt. 

so  u  bo 


(3.2.1) 
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Consequently  the  conditional  density  of  x  given  y  h(x)  is  h(x). 

When  a  pair  of  independent  uniform  (0,1)  numbers,  u  and  v,  arc 
generated  and  the  resulting  x  and  y  satisfy  the  inequality  y  <_  h (x) , 
x  is  then  a  required  random  variate  from  h(x).  Otherwise  the  pair  is 
rejected  and  a  new  pair  is  generated  and  checked.  It  is  clear  that 
the  inequality  will  be  satisfied  with  a  high  probability  only  if  the 
value  c  •  (Ci-^2)  close  to  one.  If  not,  the  procedure  will  produce 
a  large  proportion  of  inadmissable  (u,v)  pairs  and  reduce  the  efficiency 
of  the  scheme.  Even  when  this  proportion  is  small,  the  time  required 
to  check  y  <_h(x)  is  generally  rather  long  compared  with  that  of  the 
table  look-up  procedures.  Consequently,  the  acceptance-rejection 
method  is  best  used  for  generating  those  infrequent  values  from  the 
"triangular"  and  tail  regions  described  in  section  3.1. 

3.3  Generation  of  Numbers  with  a  Triangular  Density 
Let  u  and  v  be  independent  uniform  (0,1)  random  numbers.  Then 
T  *=  min  (u,v)  has  the  distribution  function 

G(t)  =  1  -  P [T>t]  =  1  -  P[u>t,  v>t]  =  1  -  (1-t)2,  0<t<l,  (3.3,1) 

and  the  density  function 

g(t)  =  2(l-t),  0<t<l.  (3.3.2) 

See  figure  3.3a,  A  linear  transformation  T’  =  a  +  bT  produces  a  random 
variable  with  density 


g(t’)  =  §  ‘  p.  (f  -  a). 


a<t<a+b 


(3.3. 3) 
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3.4  The  Distribution  of  the  Minimum  of  a  Random 
Number  of  Uniform  (0,1)  Random  Variates:  Some 
Special  Results  for  the  Exponential  Distribution 

Let  x  =  c  •  min(uj,  uj ,  . ..,  u^)  where  the  u^  are  random  indepen¬ 
dent  uniform  (0,1)  variates,  c  is  a  positive  constant,  and  N  is  a  dis¬ 


crete  valued  random  variable  with  probability  function  P  (N=n)  -  q(n) 
for  n  =  1,2,  ....  Then  the  distribution  function  of  x  can  be  expressed 


as  follows: 


F(x)  =  E  PJX  <_  x  |  N  =  n]  •  PIN  =  nj 
n=l 


=  I  [1  -  (1  -  x/c)n]  •  qfn) 
n=l 

F(x)  =  3  -  l  (1  -  x/c)n  •  q(n)  0  <  x  <  c  (3.4.1) 
n=l 

Now  consider  the  special  case 

q(n)  =  cn/[n!  (eC-l)j  n  =  1,2,  - 

The  distribution  of  x  becomes 


F(x)  =.l  -  •?  (ox)n/n!  (eC-l)  =  1  -  (eC"X-l)/(eC-i) 
n=l 


F(x)  -  (l-e~X)/(l-c'C)  0  <  x  <  c  (3.4.3) 

Thus  x  has  the  exponential  distribution  truncated  un  the  right  at  the 


value  c. 
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If,  instead  of  the  probability  function  given  in  (3.4.2),  N  has 
the  distribution  specified  by 

g(n)  =  An/n! (eA-l-A)  n  =  2,3,  ...,  (3.4.4) 

where  A  is  a  constant  as  defined  in  section  3.1,  then  the  distribution 
function  of  x  becomes 

F (x)  =  1  -  (eA-X-l-A+x)/(eA-l-A),  0<x<A,  (3.4.5) 

with  density 

f(x)  »  (e~X-e~A)/[l  -  e~A(l+A)]  ,  0<x<A.  (3.4.6) 

This  is  the  distribution  function  of  a  random  variable  from  the  "tooth" 
of  the  exponential  distribution  between  0  and  A.  See  figures  3.4(a)  and 
3.4(b). 


FIGURE  3.4(a) 


FIGURE  3.4(h) 
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Figure  3.4(c)  shows  the  density  of  a  random  variable  from  any  other  such 
"tooth"  defined  by  k  •  A  <  z  <  (k+1)  •  A. 


FIGURE  3.4(c) 


1't  is  obvious  that  all  such  densities  are  identical,  a  unique  character¬ 
istic  of  the  exponential  distribution.  Consequently,  a  random  variable 
s  t 

from  the  (k+1)'  "tooth"  is  always  produced  by  taking  z  =  k  •  A  +  x 
where  the  value  k  *  A  is  chosen  by  the  technique  described  in  section 
3.1  and  x  is  generated  by  taking  A  times  the  minimum  of  N  uniform  (0,1) 
random  numbers  where  N  has  the  distribution  (3.4.4).  Its  expected  value 
when  A  =  0.1.,,.  is  2.02  or  in  general 


E [N]  =  A  •  (eA-l)/(eA-l-A)’ 


(3.4.7) 


The  same  unique  "no  memory"  characteristic  of  the  exponential  dis¬ 
tribution  used  above  is  also  used  for  obtaining  random  variates  from 
g3,  the  tail  of  the  exponential  density  where  t  >  c. 

Let  Y  =  W  +  X  where  IV  is  discrete  valued  with  probability  function 
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PJW  =  w  s  c  *  k]  =  Ae  W  ^  Ae'c*k,  k  =  1,2,3,  ...,  (3.4.8) 

A  =  (1  -  e  }/e  and  x  is  independent  of  w  and  has  the  distribution  func 
tion  (3.4.3).  The  distribution  function  of  Y  can  be  written  as 

G(y)  *  PIY  <.  y  =  w  +  jc]  =  PlY  <_  w  *  c]  +  P[W  =  w,  X  <_  x]  (3.4.9) 

Since  the  distribution  function  of  W  is 

k  .  , 

P[W  <  w  =  c  •  kj  =  A  I  e'rc  =  1  -•  e~C  k 
j=l 

the  expression  for  G(y)  becomes 

G(y)  =  1  -  e“ (w*c^+  Ae~W(l-e~X)/(l-e‘c)  =  1  -  eC'y  ,  (3.4.11) 

y  >  c. 

This  is  the  required  distribution  function  for  variates  from  the  density 
e~*  truncated  on  the  left  at  y  =  c.  Thus  values  of  t  from  in  the 
exponential  case  are  generated  as  the  sum  of  a  discrete  variate  taking 
values  c,  2c,  3c,  ...  and  a  continuous  variate  from  the  truncated  expo¬ 
nential  density  on  (0,c). 

3.5  Generating  Numbers  from  a  ’Nearly  Triangular’ 

Density:  Application  to  the  Half  Normal  Density 

The  results  of  this  section  are  applied  to  the  half  normal  density 

i  ^ 

f(t)  =  m*  e'^2t 


0  <  t  s  ™ 


(3.5.1) 
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for  the  generation  of  random  norm.::  deviates  to  be  described  in  chapter 
4.  Consider  the  concave  and  conv<  'riangular  densities  (members  of  gj) 
depicted  in  figures  3.5(a)  and  3.5(b),' 


FIGURE  3.5(a) 


FIGURE  3.5(b) 
t 

s 

The  parallel  chords  and  tangents  enclosing  f(t)  are  determined  by  con¬ 
struction  as  indicated  in  the  figures.  In  both  cases,  the  inner  right 

triangle,  represents  most  of  the  area  under  the  density  F(t).  This  will 
.  *  •  •  • 

be  true  in  general  as  long  as  the  ratio  a/b  is  close  to  one.  A  random 
variate  from  the  density  represented  by  this  triangle  is  generated  as 
described  in  section  3.3.  Infrequently,  a  random  variate  from  the  shaded 
[  .  area  must  be  generated  by  the  acceptance-rejection  technique  described 

!  in  section  3.2.  Marsaglia  ct  al .  [12]  have  combined  these  two  procedures 

t  [ 


as  follows. 
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1)  Generate  independent  uniform  (0,1)  random  variables  u  and  v. 

2)  If  max(u,v)  £  a/b,  put  t  =  s  +  c  •  min(u,v) 

3)  If  not,  test  b|u-v|  fCs  +  c  •  rain(u,v)). 

If  yes,  put  t  =  s  +  c  •  min(u,v). 

If  no,  go  to  step  1  and  try  again. 

In  order  to  show  that  this  procedure  produces  a  variate  t  with  the  re¬ 
quired  density,  we  conveniently  take  s  =  0  without  loss  of  generality 
and  proceed  as  follows. 

Define 

m  =  min(u,v) 

M  =  max(u,v) 

* 

x  =*  c  •  m 

y  =  b(M-m)  =  b(u-v)  (3.5.2) 

Then  the  pair  (x,y)  is  uniformly  distributed  over  the  triangle  in  figure 
3.5(c) 


FIGIJUr  3.5(c) 


To  verify  this,  consider  the  joint  distribution  function  of  X  and  Y: 

F(x,y)  =  P[X  ji.  x,  Y  y]  =  Pjro  <_  x/c,  M  -  m  £  y/bj 

This  is  just  the  area  of  the  cross  hatched  region  in  figure  3.5(d)  which 
is  seen  to  be 

F(x,y)  =  2  ■  (y/b)  •  (x/c)  =  2xy/bc  ,  0<x<c  ,  (3.5.3) 


FIGURE  3.5(d) 


Thus,  (x,y)  is  a  point  randomly  chosen  from  the  larger  triangle  in  figures 
3.5(a)  and  3.5(b).  By  the  acccptance-rej ection  technique ,  if  this  point 
falls  within  the  smaller  triangle  of  figures  3.5(a)  and  3.5(b),  the  value 
x  =  t  has  the  triangular  density  -bx/c  +  a.  Now  the  condition  max(u,v) 

<_ a/b  is  equivalent  to: 

bM  <  a 


b (tbm)  <  a  -  bin  =  -bx/c.  +  a 
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or 

y  2L  -bx/c  +  a  (3,5.4) 

which  means  that  (x,y)  falls  in  the  small  triangle.  If  this  condition 
is  not  satisfied,  a  second  check  is  made  to  determine  whether  the  point 
(x,y)  falls  within  the  shaded  area  between  f(t)  and  y2  in  figures  3.5(a) 
and  3.5(b).  If  so,  then  clearly 

y  =  b(M  -  m)  =  b|u  -  v|  <_  f (x)  =  f(c  *  min(u,v))  (3.S.S) 

as  required  by  the  second  test.  The  acceptance-rejection  principle  in¬ 
sures  that  every  x  value  passing  at  least  one  of  these  two  tests  will 
have  the  density  f(t).  The  first  test  is  rapidly  made;  the  second  test 
involves  evaluation  of  f(x)  and  is  time  consuming.  The  first  test  will 
be  passed  with  a  probability 

P  [M  <_  a/b]  =  (a/b)2 

as  given  by  the  distribution  function  of  M  =  max(u.v).  It  is  therefore 
desirable  to  have  (a/b)2  close  to  1  to  achieve  a  short  average  execution 
time.  This  will  be  the  case  if  f(t)  is  "nearly  triangular." 

For  the  half  normal  density  the  a^/b^  ratios  are  obtained  as  follows. 

Case  (a)  "Concave  Triangles"  (Sec  figure  3.5(a)). 

Let  t .  =  i  •  c  =  s 

i 

ti+1  =  (i+l)c  =  s  +  c,  i  =  0,1,2,  ...  (3.5.6) 
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The  equation  of  the  chord  and  the  tangent  are  respectively: 

yl  =  "  bif-  +_^i 

y2  =  -  bi(x-s)/c  +  a ^ 

The  value  of  b.  is  clearly 

i 

b.  =  f(ti)  -  f(t.+1)  =  f(s)  -  f(s-+c) 

The  slope  of  both  lines,  -  b^/c,  must  equal  f'(X^).  Thus 

-b./c  =  -  /2/rr  x.e  ^xi  =  -  X.f(x.) 

1  l  i^i. 

This  equation  must  be  solved  iteratively  for  >L,  the  abscissa  of  the 
point  of  tangency.  Finally,  this  value  is  substituted  into  the  tangent's 
equation  to  give 

f(X.)  =  f(t.+1)  =  (-  b./c)(X.-s)  +  a..  (3.5.10) 

Substituting  the  expression  for  b  /c  and  solving  for  a.  gives 

i  •  1 

a.  =  fCX^Jl  +  Xi(X.-s))  -  f(s+c).  (3.5.11) 

Consequently  fdr  "concave  triangles" 

ai/bi  =  [f(*.)fl  +  X.(Xi-s)}  -  f(s+c)]/[f(s)  -  f(s+c)]  (3.5.12) 

i  • 

Case  (b)  "Convex  Triangles"  (See  figure  3.5(b)) 

In  tli i s  case  the  point  of  tangency  occurs  at  the  end  of  the  inter¬ 
val  where  x.  =  t.  ,  ~  s  +  c.  The  slope  of  both  lines  is  given  bv 
ii+I 

-  b./c  =  f’(s+c) 


(3.5.7) 


(3.5.8) 


(3.5.9) 
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or  b  =  -  cf  (s  +  c) .  (3.5.13) 

The  value  of  a^  is  simply 

f(t.)  -  f(t.  1  =  f(s)  -  f (s-i-c). 
i  i  +  l 

Consequently,  for' "convex  triangles" 

a./b.  -  [f(s)  -  f(s+c)]/I-  cf'(s+c)J  (3.5.14) 

3.6  Generation  of  Random  Variates  from 
the  Tail  of  the  Half  Normal  Distribution 

The  procedure  suggested  by  Marsaglia  et  al.  [12]  for  the  generation 

of  random  variates  from  the  tail  of  the  half  normal  distribution  is  based 

on  the  rejection  principle  described  in  section  3.2.  The  procedure  is  to 

generate  a  pair  of  independent  half  normal  variates  xj  and  x2  such  that 

the  point  (xj,  x.;0  lies  outside  the  quarter  circle  (see  figure  3.6(a)). 


* 
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This  is  done  by  generating  pairs  of  uniform  (0,1)  random  numbers 
and  u2  until  0  <  u2  +  u2  <_  1  and  setting 

'  =  /a'5  -  21n(u2  +  up  *  s4i^/(u2  +  u2) 

x2  =  A2  -  2 In  (Uj  +  u£)  -  /u y  (u2  +  u^} 

To  show  that  the  point  (xj ,x2)  lies  outside  the  quarter-circle,  let 
w  =  -  ln(u^  +  u|)  where  0  <  u2  +  u2  <_  1.  Then  take  the  sum 

+  xij;  =  a2  +  2w  where  0  <_  w  <  »  so  that  x2  +  x2  >_  a2.  If  the  point 

(xj.xj),  generated  in  this  fashion,  lies  within  the  square  (figure  3.6(a)), 
it  is  rejected  and  a  new  pair  of  uniform  (0,1)  random  numbers  is  generated 
and  tested.  Otherwise  the  variable,  xi  or  x2,  whose  value  exceeds  a  is 
taken  as  the  required  random  number. 

In  order  to  show  that  this  procedure  produces  the  desired  result, 
let  u  and  u  be  a  pair  of  uniform  (0,1)  random  numbers  conditioned  by 
0  <  u2  +  u2  <_  1.  Then  the  joint  density  of  Ui  and  u2  is 


f (Uj ,u2) 

=  4/tt,  0  <  U]  <  1,  0  < 

u2  <  1, 

0  <  u^  +  u2  <_  1. 

Now  define: 

(3.6.3) 

X  =  u^  +  u2,  0  <  X  < 

1 ,  and 

(3.6.4) 

Y  =  u  2  /  u  j ,  0  <  Y  < 

CO 

(3.6.5) 

Then  the  probability  distribution  of  X  is 

r(x)  =  P(X  <_  x)  =  (TT  •  x / -1 ) /  (n  ■  12/-1)  = 


(3.6.1) 

(3.6.2) 


x 


(3.6.6) 
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and  the  density  of  x  is 

f(x)  =  1,  0  <  x  <  L,.  _  (3.6.7) 


Consequently  X  is  uniform  on  (01).  The  joint  density  of  X  and  Y  is: 

g(xjy)  =  f(u,v)  .  I J  I  -  (4/u)  [1/2  (1+y2)  j  -  2/t.  ( 1+y2)  =  g^x)  •  g^y). 

(3.6.8) 

lienee  X  and  Y  a i' c  stocliuSt-icall^  i r.ucpcndsiit..  Substilut ion  o"£"  £3,6,4) 
and  (3.6.5)  into  (3.6.1)  and  (3.6,2)  yields: 


Xj  =  '^r~~2YnX  •  /l/(l+Y2) 
x2  =  /a-  -  21rTx  •  AVd+Y^). 
The  inverse  transformation  is; 

x .  ■ 

V  »  X2/X, 

The  joint  density  of  Xj  and  is  given  by: 
f(x1,x2)  =  g(x,y)  •  | J | 

I  1  , 

-  [2/ ir(l+y2)  ] (1+y2)  •  e_!5(xi+x2'a2) 

=  (i£77  e'1^1)  (/>“  e"*iX2)  caZ/2 


(3.6.9) 

(3.6.10) 


(3.6.12) 


(3.6. 12) 


I  (3.6.13) 


Clearly  x(  and  x  are  independent  random  variables  from  the  tail  of  the 
half  normal  distribution  (Xj,.\  >  a). 


CHAPTER  4 

GENERATION  OP  PSEUDO-RANDOM  NUMBERS  FROM  THE 
EXPONENTIAL,  GAMMA,  AND  NORMAL  DISTRIBUTIONS 

Before  presenting  the  procedures  for  the  generation  of  random  numbers 

from  the  exponential,  gamma,  and  normal  distributions,  a  simple  example 

is  given  demonstrating  the  use  of  the  techniques  discussed  in  Chapter  3 

for  obtaining  random  variates  from  any  continuous  distribution  defined  on 

a  finite  interval.  In  the  interest  of  clarity  all  numbers  are  expressed 

to  the  base  10.  Consider  the  distribution  defined  by 

t 

!■  (>’)  =  (y/2)  (3-y2 )  0  <  y  <  1  (4.0.1) 

with  density 

f (y)  =  (3/2)(l-y2)  0  <  y  <  1.  (4.0.2) 

Using  the  composition  technique  presented  in  3.1,  f(y)  is  represented  as 
a  mixture  of  2  densities  (since  there  is  no  tail): 

f(y)  =  PjgjCy)  +  p2g2(y).  (4.0.3) 

i 

ms  before,  g  is  a  series  of  rectangles  and  g2  represent:  the  nearly  tri¬ 
angular  regions  between  g  and  f.  See  figures  4.0(a),  (b) ,  and  (c)  . 
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As  described  in  3.1,  define  discrete  variables 

=  (i-1)  *  A  for  i=l,2,...,10  and  in  this  example  6=0.1^^^.  (4.0,4) 
Now  assign  each  a  probability  as  follows: 

Pi  =  A  ‘  fCxi+l)‘  (4.0.5) 

Now  define  a  value  Tb  as  given  in  (3.1.7),  the  area  of  the  triangular 
region  above  the  i^'  rectangle: 

Ti  *  F(xi+13  ‘  F{xi)  "  V  (4-°-6) 

These  values  are  recorded  in  table  4.0.1.  At  this  point  it  is  obvious 
that  the  values  of  pj  and  p2,  defined  in  (4.0.3),  are  pj  =  0.9225  and 
p?  s  0.0775.  Next,  the  four  sets  of  the  discrete  values  are  stored, 
based  on  the  high  order  of  digits  of  the  PJs. 
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TABLE  4.0.1 
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Since  the  P^'s  have  only  4  digits,  application  of  this  technique  will 
account  for  all  of  g  and  there  is  no  rectangular  residue  in  this  simple 
example.  Since  the  composition  of  these  4  sets  is  constant  (the  distri¬ 
bution  is  parameter  free),  the  sets  may  be  overlapped  and  considerable 
storage  space  saved.  Looking  at  table  4.0.2,  set  3  occupies  locations 
2  through  6,  set  2  occupies  locations  4  through  40,  set  3  occupies  locations 
18  through  67,  and  set  4  occupies  locations  46  through  70.  Locations  0 
through  9  are  also  used  for  the  generation  of  variates  from  g2  ,  as  de¬ 
scribed  subsequently. 


TABLP  4.0.2 
A-  *  /\bLL 


Con 

I  oc 

Con 

hoc 

Con  I 

1.0  c 

Con 

1  I.oc 

Con 

Loc 

Con 

Loc 

Con 

0.9 

10 

0,6 

20 

0.7 

30 

O.S 

40 

0 

50 

0.3 

60 

0.2 

0.8 

1.1 

0.6 

21 

0.7 

31 

0.3 

41 

0.8 

51 

0.6 

61 

0.2 

0.3 

12 

0.6 

22 

[  0.7 

32 

1  0.3 

: 

42 

0.3 

52 

0.6 

62 

0.2 

0.1 

13 

0.6 

23 

0.7 

33 

0.2 

i 

43 

0.3 

53 

0,6 

63 

0 

0.4 

14 

0.6 

24 

0.6 

34 

0. 1 

44 

0.3 

54 

0.6 

64 

0 

0.2 

15 

0.5 

25 

0.5 

35 

0. 1 

4  5 

0.3 

55 

0.6 

65 

0 

0 

16 

0.5 

26 

O.S 

36 

0.1 

46 

O.S 

So 

0.4 

66  j 

0 

0,7 

17 

0.2 

27 

0.5 

37 

0.1 

47 

0.5 

57 

0.4 

6  7 

i 

0 

0.6 

18 

O.S 

28 

0.5 

!  38 

0 

48 

0.8 

SS 

0.2 

i 

6S 

.0.4 

0.5 

19 

0.8 

29 

0.  5 

39 

0 

49 

O.S 

59 

0.2 

69 

0.4 

1 

1 

70 

1 

0.4 

1 
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As  presented  in  Chapter  2,  define: 


Qo  =  °j 

-r  10 

Q  =  10  •  26. 

i-1  ri 

(4.0.7) 

and 

no  =  o. 

r  10 

11  =2  2  5..  for  r=  1,2 ,3,4. 

r  ....  ji 
j - 1  r=l 

4 

(4.0.8) 

Using 

the  above  definitions 

define : 

s0  -  o. 

k  k 

Sk  =  10K  2  Qr 

r=l 

(4.0.9) 

and 

Nk  “  nk-l  *  10  *  Sk-l  for  1{s,1»2,3,4.  (4.0.10) 

The  N^'s  are  identical  to  the  quantity 

"s-l  *  6’  •  Pr  ' 
r=0 

in  (2.0.3).  Because  of  overlapping  between  the  four  sets,  the  N^'s 
must  be  adjusted.  Their  values  and  the  values  of  the  S^'s  are  recorded 


in  tabic  4.0.3. 
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The  procedure  for  generating  random  numbers  from  gj  is: 

1)  Generate  a  uniform  (0,1)  random  number  u  =  0.djd2.... 

2)  If  d]  <  S(i) ,  set  V  =  A(dj  +  NCI"))  +  0.1  •  (.d^dg...); 
otherwise  go  to  (3). 

3)  If  djdj  <  S (2) ,  set  Y  =  A{did2  +  N(2)}  +  0.1  •  (,d5d6...); 
otherwise  go  to  (4). 

4)  If  d!d2d3  <  S  (3) ,  set  Y  =  Atd^da  +  N(3)}  +  0.1  *  (.d5d6...); 
otherwise  go  ta  (5). 

5)  If  d1d2d3dli  <  S (4) ,  set  Y  =  A(d!d2d3d4  +  N(4))  +  0.1  •  (.d5d6...); 
otherwise  generate  a  random  variate  from  g2. 

In  order  to  generate  random  variates  from  g2»  define: 

C(l)  =  S4  +  T  £  A ( 1 ) ]  and  (4.0.11) 

C(k+1)  =  C(k)  +  T[A(k)]  for  k=l,2,...9.  (4.0.12) 

This  is  identical  to  the  definition  given  in  (3.1.13),  except  that  in 

this  example  there  is  no  rectangular  residue.  The  notation  T[A(k)]  denotes 

the  area  of  the  triangular  region  (the  T^'s  in  table  4.0.1)  corresponding 

til 

to  the  discrete  value  stored  in  the  x  location  of  the  A-table  (Table 
4.0,2).  The  values  of  the  C(k)'s  are  given  in  table  4.0.4.  If  u  is  a 
uniform  (0,1)  random  number  such  that  u  >_  S4,  the  proper  triangle  of  g2 
may  he  chosen  simply  by  testing: 

u  <  C(k)  until  the  condition  is  satisfied  (note  that  1.0) 

The  discrete  value  occupying  location  k  of  the  A-table  denotes  the  correct 
triangle. 

Having  selected  the  proper  triangle  of  g2,  a  random  variate  having 
the  density  ol  that  triangle  is  to  be  generated.  This  is  accomplished 
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through  the  use  of  the  acceptance- rejection  principle  discussed  in  3.2 
and  the  technique  given  in  3.5.  Consider  a  particular  triangle  or  tooth 
from  g2  (figure  4.0(d)).  The  tangent  and  the  chord  are  constructed  as 
indicated  in  3.5 


i 


The  inner  right  triangle  encloses  most  of  the  area  under  f(x).  A  random 
variate  from  the  density  of  this  triangle  is  generated  as  described  in 
3.3.  Occasionally,  a  random  variate  from  the  shaded  region  must  be  gener¬ 
ated  by  the  acceptance-rejection  technique.  The  values  of  a^,  b^,  and 
d.  are  obtained  as  follows: 

l 

'  -'(a./0.1)  =  f'C-V  =  3X.  (4.0.13) 

where  is  the  abscissa  of  the  point  of  tangcncy  and 


a.  =  f (x. )  -  f (x.  ) , 

x  x  x+1 


(4.0. 14) 


Substituting  (4.0.14)  into  (4.0.13)  and  solving  for  X^  yields: 

X.  =  [f(x.)  -  f(x.+1)]/0.3  (4.0.15) 

The  ordinate  of  the  point  of  tangency  is: 

y.  =  f(X.)  -  f(x.+1).  (4.0.16) 

At  x  =  x.,  >']  =  yi  and 

bi  =  yi  +  -  X.)  (4.0.17) 

The  value  of  d^  may  be  obtained  by  noting  that: 

(a^/0.1)  =  (b./d^j  and  hence  d^  =  (0.1  ♦  b^/a^.  ,  (4.0.18) 

The  values  of  a^/b.,  and  are  recorded  in  table  4.0.4  in  the 
R-table,  D-table,  and  B-table  respectively. 


C- Table 


Loc 


R- Table 


R-  Table 


B-Tab  1  o 


Con 

Loc 

Con 

Loc 

Con 

Loc 

Con 

.9370 

0 

.9870 

0 

.1013 

0 

.28875 

.9500 

1 

.9855 

1 

.  1015 

1 

.25875 

.9555 

2 

.9655 

2 

.  1036 

2 

.10875 

.9580 

•  3 

.9231 

3  4 

.  1083 

3 

.04875 

.9650 

4 

.9730 

4 

.1028 

4 

.13875 

.9690 

5 

.9524 

5 

.1050 

5 

.07875 

.9700 

6 

.8000 

6 

.1250 

6 

.01875 

.9815 

7 

.9836 

7 

.1017 

7 

!  22875 

.9915 

8 

.9811 

8 

.1019 

8 

.19875 

1.0000 

9 

.9778 

9 

.1023 

9 

.  16875 

Table  4.0.4 
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The  procedure  i'or  obtaining  random  variates  from  g2  is: 

1)  Let  u  be  a  uniform  (0,1)  random  number  such  that  u  >^4. 

2)  Test  u  <  C (k)  for  k  =  0,1,  , 9  until  the  inequality  is 
satisfied. 

3)  Generate  independent  uniform  (0,1)  random  numbers  uj  and  u2. 

4)  If  iikix(u1,u2)  JL  R(k),  then  set  Y  =  A(k)  +  D(k)  *  min(u1}u  ); 
otherwise  go  to  (5) . 

5)  If  B(k)  •  |uj  -  u2|  <_  f[A(k)  +  D(k)  *  min(u1>u2)], 

set  Y  =  A(k)  +  D(k)  •  min(u1(u2);  otherwise  go  to  step  (3) 
and  try  again. 

A  schematic  of  the  entire  procedure  for  generating  random  variates 
from  the  density  f(y),  given  in  (4.0.2),  is  presented  in  figure  4.0(e). 

4.1  Generation  of  Exponentially 
Distributed  Random  Numbers 

Applying  the  composition  principle  presented  in  3.1,  the  exponential 
density  function, 

f(t)  e-t,  0  <_  t  <_  (4.1.1) 

*  ■  /  . 

is  represented  as  a  mixture  of  three  densities 

f(t)  =  P^Ct)  +  P2g2(t)  +  p3g3(t).  ;  (4.1.2) 

As  before,  g  is  a  scries  of  rectangles  defined  for  t  <_  4,  g2  represents 
the  toothlike  region1'  between  g  and  f  defined  for  t  £  4,  and  g3  is 
the  tail  of  f  for  t  >_  4.  The  probabilities  p  ,  p2,  and  p3  arc  rcspcc- 


S3 


tively  the  areas  represented  by  g,,  g  ,  and  g 3 . 

As  described  in  3.2  a  random  variable  Y  having  density  gj  is  generated 
as  follows: 

Define 


M  =  4/A,  where  A  =  10  *  for  a  decimal  machine  or  A 
binary  machine;  and  discrete  variates 
=  (i-1)  •  A  for  i  =  1,2,  ...,  M. 

Each  is  assigned  a  probability 


P. 

l 


A  •  f  (t .  .  )  =  A 
v  l  +  l 


6  .6  .6  .  , 

li  n  3i 


16  ^  for  a 

(4.1.3) 

(4.1.4) 


(4.1.5) 


Now,  using  the  high  order  4  octal  digits  (binary  machine)  or  the  high 
order  3  decimal  digits  (decimal  machine),  store  the  t^'s  according  to  the 
technique  described  in  Chapter  2.  As  in  the  case  of  the  example  given  in 
the  beginning  of  this  chapter,  the  sets  of  discrete  variates  are  stored 
in  the  manner  that  permits  maximum  overlap.  See  the  D-table  in  the  pro¬ 
gram  listing  of  GEN1  in  the  appendix;  set  1  occupies  locations  49A  thru 
4BE,  s  2  occupies  locations  40C  through  4AC,  and  set  3  occupies  locations 
37D  through  4S3,  Locations  374  through  3B3  are  used  for  the  generation 
of  variates  from  g2  and  the  residual  of  as  described  in  3.1. 

At  this  point  define: 

M 

%  =  Q,  •  E  5  (4.1.0) 

i-1 


r  M 

H0  =  0,  n  -  E  E  6..  for  r  =  1,2, 3, 4  (4,1.7) 

j=l  i-1 

for  binary  machines  where  the  6's  arc  octal  digits  and  C  =  S  or  r  =  1,2,3 
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for  decimal  machines  where  the  6's  are  decimal  digits  and  6  -  10.  Also 
define : 


S0  -  o,  sk 


(4.1.8) 


and 

Nk  =  nkl  -  &  •  Sw_j  for  k  =  1, 2,3,4  (4.1.9) 

in  the  binary  case  (8  =  8)  or  k  =  1,2,3  in  the  decimal  case  (8  =  10). 

The  S-values  and  the  N-values  are  recorded  in  the  program  listing  of  GEN1 
in  the  appendix.  Note  that  the  N-values  have  been  adjusted  to  allow  for 
overlap  between  the  sets  of  the  D-table. 

\ 

The  procedure  for  the  generation  of  random  variates  from  this  trun¬ 
cated  portion  of  gj  is:  (binary  case) 


1)  Generate  a  uniform  (0,1)  random  number  u  *  .d^d  .... 

2)  If  d1dJ>  <  S2>  set  Y  =  D(djd2  +  N2)  +  A  •  (.d^...); 
otherwise  go  to  (3). 

3)  If  d1d2d3  <  S  3 ,  set  Y  =  D^d1d2ci3  +  NV  +  A 
otherwise  go  to  (4) . 

4)  If  d  d  dd,  <  S,  ,  set  Y  =  D(d,d,d,d.  +  N,  }  +  A  •  (.d  d  ...); 
otherwise  generate  a  random  variate  from  the  residual  of  g(,  g2, 
or  g3. 

Note  that  in  the  above  procedure  the  testing  begins  with  S2.  This  is 
because  5  ^  =  0  for  all  i  since  0  <  f(x)  <  1  for  0  <  x  <  4. 

Random  variates  from  the  residual  of  gj  and  from  are  generated 
by  first  selecting  one  of  the  rectangles  of  g  with  appropriate  probabil¬ 
ity  and  then  testing,  to  determine  whether  the  required  variate  should 


come  from  the  tooth  or  from  the  residual  corresponding  to  that  rectangle. 
This  is  accompl i sited  as  follows: 

In  M  consecutive  locations  store  each  of  the  t.’s  once  (locations  374 

i 

through  3B3  of  the  D-tablc  ir«  the  listing  of  GFN’l  in  the  appendix).  For 
each  define  a  pair  of  values  as  given  in  (3.1.12),  (3.1.13),  and  (3.1.14) 


C(l)  =  s,  ♦  R[D(1)J, 

(4.1.10) 

E(k)  =  C(k)  +  T[D(k)]  and 

(4.1.13) 

C(kU)  =  C(k)  +  R[D(k  +  l)  ] . 

(4.1.12) 

In  the  exponential  case 

T[D(k) ]  =  e'D(k;i  -  (1  +  A)  .  e’ (DCk)+A:) 

(4.1.13) 

and 

R[D(k)J  =  P[D(k)]  -  .<5[l,D(k)j6[2,D(k)]6[3,D(k)]6[4,D(k)),  (4.1.14) 

where  P[D(k)]  is  as  defined  in  (4.1.5).  The  values  T[D(k'J]  and  R[D(k)] 
denote  respectively  the  areas  of  the  triangular  region  and  rectangular 
residue  lost  upon  truncation  of  PfD(k)]  corresponding  to  the  t-value 
occupying  location  D(k).  The  values  of  the  B(k)'s  and  the  C(k)’s  are 
also  recorded  in  the  listing  of  Gh.N'l  in  the  appendix.  A  random  variate 
with  a  rectangular  distribution  is  generated  as  described  in  5.1,  and  in 
the  exponential  case  a  variate  having  the  density  of  one  of  the  teeth  of 
#2  is  generated  by  the  technique  given  in  3.4.  The  procedure  for  gener¬ 
ating  random  variates  from  either  the  residual  of  gl  or  from  g2  is  sum¬ 


marized  as  iollows: 
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1)  Let  u  be  a  uniform  (0,1}  random  number  such  that 

Mu<  Pj  +  P2  =  1  -  c’'*- 

2)  Test  u  <  1S(J)  for  J  =  1,2,  M,  until  the  inequality  is 

satisfied  (note  B(M)  =  1  -  e  **)  for  same  J  =  J'.’ 

3)  !f  u  <  C(J'),  generate  a  new  uniform  (0,1)  random  number  v  and 
set  Y  =  D(-J')  +  a  •  v  (this  produces  a  number  from  gj); 
otherwise  go  to  (4). 

4)  Generate  a  new  uniform  (>'  ■  random  number  v  and  test 

n 

v  <  Q(n)  =  Z  A3/[j!(e  -  1  -  &)]  for  n  =  2,3,  ... 

j  =  2 

until  the  inequality  is  satisfied  for  n  =  n1  . 

5)  Generate  n*  independent  uniform  (0,1)'  random  numbers  and  set 
Y  =  D(J')  +  A  •  min(u1,u2,  ....  u^,).  This  produces  a  number 
from  g2. 

A  random  variate  Y  having  density  is  generated  as  a  sum 

Y  =  W  +  X  (4.1.15) 

where  X  is  exponential  on  (0,4)  and  IV  has  the  distribution 

P[W  =  w  =  4  •  k]  =  e"*4^ (e4  -  1),  k  =  1,2,  _  (4.1.16) 

♦  ’  •  . 

The  procedure  for  generating  a  random  variate  from  is  summarized  as 
follows: 

1)  Let  u  be  a  uniform  (0,1)  random  number  such  that  u  >_ pj  +  p2  = 

1  -  c'4. 

-‘i  (k+P 

2)  Test  u  <  1  -  c  '  for  k  =  1,2,  ...  until  the  inequality  is 

satisfied  for  some  k  =  k'. 


3)  Set  W  -  4  *  k1  and  gene-rare  a  new  uniform  (0,1)  random  number  v. 

- 1» 

4)  Set  u  =  (pj  +  p2)  •  v  =  (1  -  c  )  •  v  and  enter  the  procedure  for 

generating  exponentially  distributed  random  variates  on  the  inter- 

—  4 

val  (0,4)  with  u  --  (1  c  )  •  v  as  the  initial  uniform  number 
and  generate  X. 

5)  Set  Y  =  IV  +  X. 

In  [13]  Marsaglia  presents  a  schematic  and  the  necessary  constants 
for  the  generation  of  random  variates  from  the  exponential  distribution. 
Application  of  this  method  results  in  a  very  fast  generator  program.  As 
is  evident  from  the  discussion  herein  the  method  is  also  exact;  the  accur¬ 
acy  of  the  result  is  dependent  only  on  the  word  sice  of  the  computer. 

D 

4.2  Generation  of  Gamma  Distributed  Random  Variates 
In  order  to  produce  random  variates  having  the  gamma  distribution, 
the  following  result  is  used. 

If  is  a  random  variable  with  density 

-x. 

f(x.)  =  e  *,  0  <  “  for  i  =  1,2,  . n  (4.2.1) 

then  the  random  variable 

l 

n 

Y  -  lx.  (4.2.2) 

i*l  1 

has  the  ganuma  density 

g (y)  =  (l/r(n))yn  1cy,  0  <  y  <  «  (4.2.3) 

when  the  x.'s  are  stochastically  independent.  Consequently,  the  genera¬ 
tion  of  random  variates  from  the  gamma  distribution  involves  only  a  simple 


Huu  a  i  i  ^  u t.  « it  u  1  t  _  ,  ",  one...  1  ;j  1  ;  rt  I’.'jd  , 

1)  Us  ini;  the  net  hod  g'.vi-n  in  1.1,  generate  r.  exponentially  distributed 

random  numbers  x  .  ,  x, .  ....  x  . 

•  n 

n 

2)  Set  Y  =  S  x.  (if  a  scale  parameter  b  t  1  is  desired,  set 

i  =  l  1 
n 

Y  =  B  •  i  x.).  * 

i=l  1 


4.3  Generation  of  Random  Variates 
having  the  Standard  Normal  Distribution 

Using  the  composition  principle  of  3.1,  the  absolute  normal  density 

f(t)  =  /2/7  e^t2,  t  >_  0  (4.3.1) 

is  represented  as  a  mixture  of  three  densities 

f(t)  =  Pj£j(t)  +  p2g2(t)  +  p3g3(t).  (4.3.2) 

Again,  gj  is  a  series  of  rectangles  defined  for  t  <_  3,  g2  represents  the 
toothlike  regions  above  g  and  below  f  defined  for  t  <_  3,  and  g3  is  the 
tail  of  the  density  (4.3.1)  defined  for  t  >_  3. 

Random  variates  from  gj  are  generated  by  means  of  the  technique  given 
in  3.1.  Define 

M  =  3/e  where  A  =  10  *  (decimal  machine) 

or  &  =  16  1  (binary  machine),  and  (4.3.3) 

discrete  variables 

=  (i  -  1)  •  A  for  i  =  1,2,  ...,  M. 


(4.3.4) 


Ass  i  i'll  each  t  a  :’r 

l  ‘ 


P.  =  A  •  f(t  .  ,  ' 

1  s  J+] ■ 


2/~  c 


i  *  i 


f  f, 

J  l  / 1 


t*  ..'.5] 


Applying  the  storage  technique  described  in  Chapter  2,  store  the  t  's 

based  on  the  high  order  4  octal  digits  (binary  machine)  or  the  high  order 

« 

3  decimal  digits.  Note  that  as  in  the  exponential  case 


i  i 


0  for  all 


i,  and  that  the  sets  of  discrete  values  are  overlapped  (see  the  A-table 
in  the  listing  of  GEN3  in  the  appendix).  Now  define  values 

M 


and 


Jlso  define 


and 


Qo  =0,  Q  =  6  £  fi 

i=l  r 


r  M 

JI0  =  0,  n  =  £  £  5  • 

r  j-1  i=l  J1 


s0'  =0,  S  =  BK  •  £  Q 

k  r=l  1 


Ni  =  nk-l  '  6  ’  Sk-1' 


(4.3.6) 


(4.3.7) 


(4.3.8) 


'  (4.3.9) 


for  a  binary  machine  r,k  =  1,2, 3,4  and  6=8  for  a  decimal  machine, 
r,k  =  1,2,3  and  6  =  10.  The  S-values  and  N-values  are  also  recorded  in 
the  program  listing  of  G11N3  (note  that  the  N-values  have  been  adjusted  to 
allow  for  overlap  between  sets). 

The  procedure  for  the  generation  of  random  variates  from  this  trun¬ 
cated  portion  of  gj  is:  (binary  case) 

1)  Generate  a  uniform  (0,1)  random  number  u  =  . did2.... 

2)  If  d  ]  d2  <  S?,  set  Y  =  A  { d  j  cl  2  +  \’2}  +  A  *  (.d5d6...); 


otherwise  go  to  (3). 
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3)  If  <  S^,  set  Y  =  A{d1d2d3  ■*  Nj)  +  A  •  (.d&d£...)j 

otherwise  go  to  (4). 

4)  If  djd2d3d(;  <  S4.  set  Y  =  Afd^d.d^  +  N4)  +  A  •  (.dgdg...); 

otherwise  generate  a  random  variable  from  g^,  g^,  or  the 
residual  of  £,. 

Random  variates  from  the  residual  of  gj  and  from  g^  are  produced 
by  first  choosing  one  of  the  rectangles  of  g  with  appropriate  probability 
and  then  testing  to  determine  whether  the  required  variate  should  come 
from  the  rectangular  residue  of  or  from  the  toothlike  region  of  g2. 

As  presented  in  3,1,  this  is  accomplished  as  follows: 

In  M  consecutive  locations  store  each  t^  once,  and  for  each 
define  a  pair  of  values  as  ii  (3.1.12),  (3.1.13),  and  (3.1.14): 

C(J)  =  S4  ♦  R[A(l)]f  (4.3.10) 

B(10  =  C(k)  +  T[A(k)] ,  (4.3.11) 

C(k+1)  =  C (k)  +  R[A(k+l)] ,  (4.3.12) 

Note  that  in  the  listing  of  GEN3  this  portion  of  the  A-table  occupies 
locations  690  thru  6BF.  However,  for  convenience  it  is  assumed  that 
location  690  correspond  to  A(l).  R[A(k)]  and  T[A(k)]  denote  respective¬ 

ly  the  area  of  the  rectangular  remnant  and  the  area  of  the  toothlike 
region  corresponding  to  the  t-value  occupying  location  A(k).  For  the 
normal  case: 


R(A  (k)  ]  =  P[A(k)J  -  .i5[l,A(k)]5[2,A(k)]6[3,A(k)]6[4,A(k)] 


(4.3.13) 
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where  P[A(k)]  is  defined  in  ('1.3.5)  as 

P[A(k)]  =  A  •  f  (A  (k)  +  A)  =  A  •  /2h_  •  e'^(A(k)  +  &)2 

=  ,6[l,A(k)]<5(2,A(k)]. . .  ,  ,  £4.3.14) 

and 

T[A(k) ]  =  F (A (k)  +A)  -  l-(A(k))  -  P[A(k)]  (4.3.15) 

where  x  ,  2 

F (x)  =  JljH  J0  e'"J  dt  (4.3.16) 

Once  having  determined  whether  a  variate  with  a  rectangular  density  (g.) 
or  a  variate  with  a  toothlike  density  (g2)  is  required,  the  schemes 
given  in  3.1  and  in  3.5  are  used  to  produce  the  variate  as  required.  This 
procedure  is  summarized  as  follows: 

1)  Let  u  be  a  uniform  (0,1)  random  number  such  that  S4<_ u  <  px  +  p?  = 

F(5). 

2)  Test  u  <  B(k)  for  k  =  1,2,  ...,  M  until  the  inequality  is  satisfied 
for  some  k  =  k' . 

3)  If  u  <  C(k'),  then  generate  a  new  uniform  (0,1)  random  number  v 
and  put  Y  =  A(k')  +  A  •  v  (this  produces  a  variate  from  g{); 
otherwise  go  to  (4). 

4)  Generate  independent  uniform  (0,1)  random  numbers  Uj  and  u2 . 

5)  If  max(u1,u2)  <_D(k’)  ,  set  Y  =  A(k')  +  A  •  min(u  ,u  ).  This 
produces  a  variate  from  g2  as  described  in  3.5;  otherwise  go  to  (6) 

6)  Set  W  =  -  h(A  •  inin(Uj,u2)  -  A)  *  [2  •  A(k')  +  A  •  min(uJ,u2)  +  A] 


*  The  0(k)!s  arc  tabled  values  of  the  a./b.  ratios  given  in  (3.5.12) 
and  in  (3.5.14)  corresponding  to  the  u  storeA  in  A(k'). 
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and  test  1 11  j  -  u2|  <  i-(k')  •  (eW  -  1)  .  If  yes,  set 
Y  =  A(k')  +  A  *  min (uj » u  ) .  If  no  go  to  step  (4)  and  try  again. 

A  random  variate  from  the  density  r3  is  produced  by  the  technique 
presented  in  3.6.  An  absolute  normal  variable  |y|,  conditioned  by  j Y |  >  3, 
is  required.  The  procedure  is  summarized  as  follows: 

1)  Generate  pairs  of  uniform  (0,1)  random  numbers  Uj  and  U2  until 
the  condition  +  u^  <_  1  is  satisfied. 

2)  Form  pairs  and  x2  as  described  in  (3.6.1)  and  (3.6.2) 

Xj  =  •  [{9  -  21n(Uj  +  upj/tuf  +  up)'5 

x2  =  u2  •  [{9  -  21n(u*  +  u2)}/{u|  +  up]"5 

e 

3)  Test  x,  >  3.  If  yes,  set  Y  =  Xj .  If  no,  test  >  3,  If  yes,  put 
Y  =  x2.  If  no,  go  to  step  (1)  and  repeat  the  procedure. 

In  [12]  Marsaglia  presents  a  schematic  and  the  required  constants 
for  the  generation  of  random  variates  from  the  absolute  normal  distribu¬ 
tion.  This  method  results  in  random  variates  having  the  density  (4.3.1). 

In  order  to  obtain  random  variates  from  the  standard  normal  distribution. 

_ _  .1*4.2 

f(t)  =  (1//2tT)  e2  >-'x.<t<aj  ,  (4.3.17) 

1 

a  random  +  or  -  sign  must  be  attached  at  some  point  in  the  procedure.  In 
GEiN'3  this  is  done  by  generating  a  uniform  (0,1)  random  number  u  and  testing 
u  <  Jj.  If  yes,  a  -  sign  is  affixed  to  the  variate;  otherwise  a  positive 
variate  is  returned. 

The  procedures  presented  in  this  chapter  derive  their  speed  from  the 

This  expression  is  simply  ail  alternative  way  of  testing  b^|u-v| 
f (A(k 1 )  +  A  •  min(u  ,u  )). 
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fact  that  the  vast  majority  of  variates  generated  come  from  the  truncated 
portion  of  It  requires  only  the  generation  of  a  uniform  (0,1)  random 

number,  a  series  of  compare  instructions,  and  a  tabic  lookup.  SVhilc 
the  generation  of  variates  from  g?  and  may  he  somewhat  time  consuming, 
such  variates  are  required  only  rarely.*  Consequently,  that  average  time 
per  generation  is  quite  fast.  As  is  evident  these  methods  are  exact. 
Precision  is  limited  only  by  the  word  size  of  the  computer. 


CHAPTER  5 

RESULTS  OF  'I  LETS  PERFORMED  ON 
SEQUENCES  OT  PSEUDO- RAN'I«;M  NUMBERS 


S.l  Rando/n  Numbers  From  The 
Uniform  (0,1)  Distribution 

The  tests  performed  on  sequences  of  independent  uniform  (0,1)  random 
numbers  were  suggested  by  MacLaren  and  Marsagl ia  in  [9],  The  stringency 
of  these  tests  is  justified  in  that  the  procedures  for  the  generation  of 
random  variates  from  other  distributions  depend  heavily  upon  the  method 
used  to  obtain  uniform  (0,1)  random  numbers. 

The  tests  made  were  chi-square  (x2)  goodness  of  fit  tests  on  the 
distribution  of  the  random  numbers,  pairs  of  random  numbers,  triples  of 
random  numbers,  and  the  maximum  and  minimum  of  n  random  numbers.  In 
general  a  sequence  of  l  variables  Y^.Y^...,  Y^  was  calculated  from  the 
sequence  of  random  uniform  variates.  The  range  of  the  was  divided 
into  m  cel's  of  equal  probability,  p,  and  the  number  of  occurrences,  Cm, 
in  each  cell  counted.  The  x2  statistic 

in 

X2  =  Z  (0.  -  l  •  p)?/i  •  p  (5.1.1) 

i  =  l 

with  in  -  1  degrees  of  freedom  was  calculated  and  transformed  to  a  standard 
normal  deviate 

T  -  (2  •  x2y  -  (2  •  (m  -  1)  -  1)?5 


(5.1.2) 
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for  each  test.  This  form  is  valid  for  n  ?  31.  The  significance  level 
of  this  normal  deviate  was  then  computed  as 


T 

r 

0 


-Ft2 


S  =  fA  C 3 / ►//2ir )  e  *’t  dt  , 


(5.1.3) 


In  each  test  a  total  of  100,000  uniform  numbers  was  generated.  For 
convenience  the  numbers  were  generated  in  10,000  sets  of  10  uniform 
numbers  each.  Five  separate  runs  were  made  and  the  tests  performed  as  fol¬ 
lows  : 

1)  Uniformity.  The  unit  interval  was  divided  into  1000  segments. 

Each  segment  has  an  expected  value  of  100  occurrences.  The  x2  statistic 
was  computed  as 


1000 

X2  =  £  (0.  -  100)2/100 

i=l  1 


(5.1.4) 


.th 


where  Ch  is  the  number  of  observed  occurrences  in  the  i  segment. 

2)  Pairs.  Successive  pairs  of  uniform  numbers  were  taken  as  the 
coordinates  of  a  point  in  the  unit  square.  The  unit  square  was  par¬ 
titioned  into  100  cells,  each  with  an  expected  value  of  500  occurrences. 
The  x2  statistic  was  computed  as 


10  10 

X2  =  £  ?.  (0.  .  -  500)  2/500 

i=l  j=l  13 


(5.1.5) 


where  0. .  is  the  number  of  observed  occurrences  in  the  ij tn  cell. 

ij  J 

3)  Triples.  Successive  triples  of  uni  form  numbers  were  taken  as 
the  coordinates  of  a  point  in  the  unit  cube  (every  tenth  number  was  skip¬ 
ped).  The  unit  cube  was  partitioned  into  1000  cells  with  an  expected 
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value  of  30  occurrences  in  each  cell.  The  x?  statistic  was  computed  as 

10  10  10 

X?  =  £  t  5:  (0  -  30)2/30  (5.1.6) 

i=l  j=l  k~l  1JK 


where  0...  is  the  number  of  observed  occurrences 
ijk 

4)  Maximum  of  n  random  values.  If  Uj,u2,  . 
uniform  (0,1)  deviates,  then 


in  the  ijkt^  cell. 

. ,  are  independent 


W  =  max (u j ,u2,  . . . ,  un) 
should  have  the  distribution 


(5.1.7) 


P(W  <_  a)  =  F(a)  =  a.n  for  0  .<  a  <  1,  and  (5.1.8) 


F(W)  =  [nax(u1,u2. 


(5.1.9) 


should  be  uniformly  distributed  over  the  interval  (0,1).  A  total  cf 
10,000  IV's  were  generated  for  each  n,  and  F (K)  was  tested  for  uniformity 
using  the  x2  goodness  of  fit  test  for  100  equal  subintervals.  The  yz 
statistic  was  computed  as 

100 

X2  =  T.  (0.  -  100)2/100  (5.1.10) 

i  =  l  1 

where  0.  is  the  number  of  observed  occurrences  in  the  subinterval. 

In  this  test  n  takes  the  values  2  and  5. 

5)  Minimum  of  n  random  values.  This  test  is  the  same  as  that  for 
the  maximum  of  n  except 


W  =  min(Uj,u2,  ...,  u  )  and 


(5.1  .11) 


6S 


and 

F (W)  -  1  -  (1  -  W)11  (5.1.12) 

which  should  be  uniformly  distributed  over  the  unit  interval.  The 
.10,000  W's  were  segmented  into  100  equal  subintervals  and  the  test 
computed  as  in  (S.1.10).  In  this  tcst«n  takes  the  values  3  and  10. 

The  results  of  these  tests  are  reported  in  Table  5.1.1.  The  results 
compare  favorably  with  those  presented  in  [9]  and  justify  the  use  of  this 
particular  uniform  (0,1)  random  number  generator. 
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TABLE  5.1.1 


RUN  01 


TEST 

CHI-SQUARE 

D.F. 

SIGNIFICANCE 

UNIFORMITY 

975.400 

999 

0. 30186688 

PAIRS 

89.6720 

99 

0. 2S987495 

TRIPLES 

1034.7333 

999 

0.78918171 

MAX  OF  2 

97. 1800 

99 

0.46241028 

MAX  OF  5 

80.7800 

99 

0.09257838 

MIN  OF  3 

109. 5200 

99 

0.77766504 

MIN  OF  10 

75.0400 

99 

0.03713434 

RUN  02 


TEST 

CHI-SQUARE 

D.F. 

SIGNIFICANCE 

UNIFORMITY 

1047. 1600 

999 

0.85902187 

PAIRS 

92.8440 

99 

0.34129536 

TRIPLES 

1022.4667 

999 

0.70302930 

MAX  OF  2 

104.4800 

99 

0.66267689 

MAX  OF  5 

96.9800 

99 

0.45671365 

MIN  OF  3 

102.1000 

99 

0.60032472 

MIN  OF  10 

97.7400 

99 

0.47836695 
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TABLE  5,1,  ! 

RUN  HZ 

C,,I- SQUARE 


.TEST 

UNIFORMITY 
pairs 
triples 
Max  of  2 
max  of  s 
min  Of  3 
min  of  io 


1006. 3400 
93.8720 
1052 . 0000 
93.  5400 
104.5000 
102.  7400 
^07, 0600 


(continued) 

D.F. 

999 

99 

999 

99 

99 

99 

99 


significance 

0.56949803 
36930366 
0-88134258 
0. 36019046 
0-66318213 
0-01749 364 

c 

72480337 


TEST 

uniformity 

PAINS 

triples 

MAX  OF  2 
MAX  OF  5 
MIN  OF  3 

min  of  io 


RUN  # 4 


CH I -SQUARE 

D.P. 

070.2600 

999 

103. 9040 

99 

10 84.4667 

999 

83.2400 

99 

96.4200 

99 

76. 7000 

99 

116. 5000 

99 


significance 

0-26223433 
0.64798213 
0.97021768 
12861638 
0. 44078290 
0.04945168 

0. 89040197 
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TABLE  5.1.1  (continued) 


RUN  #5 


TEST 

CHI-SQUARE 

D.F. 

SIGNIFICANCE 

UNIFORMITY 

1052 . 4000 

• 

999 

0.88303445 

PAIRS 

82.5080 

99 

0.11705985 

TRIPLES 

1013.5333 

999 

0.63124447 

MAX  OF  2 

92.6200 

99 

0.33528034 

MAX  OF  5 

98.8400 

99 

0.50965471 

MIN  OF  3 

140.6600 

99 

0.99687357 

MIN  OF  10 

136.3000 

99 

0,99333696 
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5.2  Random  Numbers  Having 
A  Discrete  Distribution 

Sequences  of  random  numbers  from  the  binomial  density 

f(x)  =  (*)  *  pX  •  (1  -  p)N'X,  0  <  p  <  1,  x=0 , 1 ,,,  .  ,N  (5.2.1) 

were  generated  for  a  variety  of  parameter  combinations.  The  range  of 
the  variable  was  divided  into  k  intervals  of  probability  q^,  such  that 
the  expected  value  of  each  interval,  m  •  q^,  was  at  least  10.  The  x2 
statistic  was  computed  as 

k 

X2  =  l  (0.  -  m  •  q  )2/m  •  q  (5.2.2) 

i=l  ' 


tH 

where  0^  is  the  number  of  observed  occurrences  in  the  i  interval,  and 
m  is  the  size  of  the  sample  generated.  The  significance  level  was  com¬ 
puted  as 


XZ/2 

S  =  (1/T (v/2) )  • 


-t 

e 


t<v/2  -  1}dt 


(5.2.3) 


where  v  =  degrees  of  freedom,  in  this  case  v  =  k-1.  The  results  of 
these  tests  are  recorded  in  table  5.2.1.  A  sample  size  of  1000  random 
numbers  were  generated  for  each  test. 

Random  numbers  'from 'the  Poisson  density 

f(x)  =  *VX/x!  ,  0  <  A  <  »  ,  x  =  0,1,...,“  (5.2.4) 

I 

and  from  the  negative  binomial  density 

f(x)  =  (r+x_1)  ‘  Pr  *  U’P)\  r  >  0,  0  <  p  <  1,  (5.2.5) 


X  -  0  1  2  CO- 

A  ”  y  I  1  ^  J  •  *  '  I  J 
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were  tested  in  the  same  manner  as  random  numbers  from  the  binomial  density 
(5.2.1).  The  results  of  these  tests  are  recorded  in  table  5.2.2  and  5.2.3 
respectively.  .  „  — - -  — 
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TABLE  5.2.1 

RESULTS  OF  CHI-SQUARE  GOODNESS  OF  FIT  TESTS 


ON  RANDOM 

NUMBERS 

FROM  THE  BINOMIAL 

DENSITY 

p 

N 

D.F . 

CHI-SQUARE 

SIGNIFICANCE 

0. 100 

10 

4 

6.03494 

0.8034458581 

0.  100 

20 

6 

6.44102 

0 . 624360063S 

0.100 

30 

7 

2.52631 

0.0748969943 

0. 100 

40 

9 

6.33123 

0.2936355148 

0. 100 

50 

9 

12.92191 

0.8338284965 

0.200 

10 

5 

7.56431 

0.8180614798 

0.200 

20 

8 

16.34847 

0.9623445588 

0.200 

30 

10 

15.13733 

0.8728654417 

0.200 

40 

11 

9.00770 

0.3788180126 

0.200 

SO 

13 

12.00444 

0.4727203102 

0.300 

10 

7 

7.97911 

0.6655606411 

0.300 

20 

9 

10.98540 

0.7232888983 

0.300 

30 

11 

17.78713 

0.9133467551 

0.300 

40  ' 

13  . 

12.08372 

0.4792085395 

0.300 

50 

15 

21.66071 

0.8829858823 

0.400 

10 

7 

7.29406 

0.6010818180 

0.400 

20 

10 

12.18194 

0.7269356404 

0.400 

30 

12 

20.64248 

0.9441325898 

0.400 

40 

14 

10.05124 

0.2415723570 

0.400 

50 

16 

12.35961 

0.2811236434 

0.500 

10 

8 

6.86156 

0.4483592431 

0.500 

20 

10 

12.57912 

0.7518355135 

0.500 

30 

1.2 

10.41700 

0.4205727006 

0.500 

40 

14 

23.28293 

0.9441810635 

0.500 

50 

16 

9.42171 

0.1049637735 

0.600 

10 

7 

6.92477 

0.5632454466 

0.600 

20 

10 

18.27415 

0.9494S85520 

0.600 

30 

12 

10.17474 

0.3993657394 

0.600 

40 

14 

11.25580 

0.3341658613 

0.600 

50 

16 

10.03645 

0. 1352S1S684 
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TABLE 

5.2.1  (continued) 

\  p 

N 

D.F. 

CHI-SQUARE 

SIGNIFICANCE 

: 

0.700 

10 

7 

3.47126 

0.1617396643 

0.700 

20 

9 

10.70221 

0.7033259190 

0.  700 

30 

11 

12.02199 

0.6380075901 

0.700 

40 

13 

13.88283 

0.6178529373 

>  0.700 

50 

15 

18.60988 

0.7680404078 

;  0.800 

10 

5 

4.53985 

0.5254310881 

'  0.800 

20 

8 

16.74590 

0.9671335640 

0.800 

30 

10 

10.97150 

0.6402573630 

;  0.800 

40 

11 

9.02074 

0.3800224541 

i  0.800 

50 

13 

15.05620 

0.6961423824 

0.900 

10 

4 

3.31137 

0.4928666476 

0.900 

20 

6 

4.40204 

0.3775593406 

0.900 

30 

7 

6.78871 

0.5487935778 

0.900 

40 

9 

6.74965 

0.3368342656 

0.900 

50 

9 

5.04256 

0.1694175512 

RESULTS  OF 

RANDOM  1 

Cli I -SQUARE  GOODNESS  OF 

NUMBERS  FROM  THE  POISSON 

FIT  TESTS  ON 

DENSITY. 

'  LAM 

D.F 

CHI-SQUARE 

SIGNIFICANCE 

5.00 

10 

13.57002 

0.8068012578 

10.00 

15 

8.51890 

0.0986988385 

15.00 

17 

16.33853 

0.5000245564 

20.00 

20 

18.35403 

0.4359007928 

25.00 

21 

28.70494 

0.8787325544 

30.00 

23 

24.94425 

0.6468336304 

35.00 

25 

26.91360 

0.6397852204 

40.00 

27 

43.19749 

0.9750174127 

45.00 

28 

22.55633 

0.2450911949 

50.00 

30 

23,26187 

0.1956567510 

55.00 

30 

42.01434 

0.9286283430 

60.00 

31 

34.81330 

0.7087406320 

65.00 

33 

29.49011 

0.3573849816 

70.00 

34 

29.73472 

0.3231975827 

75.00 

35 

26. 18810 

0.1410044103 

80.00 

35 

26.80832 

0. 1619395204 

85.00 

36 

21.61020 

0.0278060381 

90.00 

37 

28,29852 

0.1527661009 

95.00 

38 

30.51639 

0.1992224404 

100.00 

39 

37.10S50 

0.4435994584 

TABLE  5.2.2 
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TABLE  5.2,3 

RESULTS  OF  CHI-SQUARE  GOODNESS  CF  FIT  TESTS  ON  RANDOM 
NUMBERS  FROM  THE  NEGATIVE  BINOMIAL  DENSITY. 


p 

R 

D.F. 

CHI-SQUARE 

SIGNIFICANCE 

t 

0.100 

2 

41 

36.08627 

0,3115372475 

0.100 

4 

57 

51.64328 

0.3244292796 

0.100 

6 

67 

77.44254 

0.8201664431 

0.100 

8 

72 

54.21390 

0.0583581269 

0.100 

10 

76 

58.52535 

0.0683S03608 

0.200 

2 

23 

20.00081 

0.3581367902 

0.200 

4 

33 

27.73007 

0.2730820886 

0.200 

6 

40 

43.00986  ■ 

0.6563763897 

0.200 

8 

45 

42.83166 

0. 43S7673160 

0.200 

10 

49 

38.69528 

0.1455417220 

0.  300 

2 

IS 

11.10495 

0.2548804445 

0.300 

4 

22 

25.94653 

0. 746015163S 

0.300 

6 

28  ' 

17.46899 

0.0612580512 

0.300 

8 

31 

25.10221 

0.2368731801 

0.300 

10 

34 

18.65952 

0.0151991805 

0.400 

2 

11 

12.32864 

0.6605459138 

0.400 

4 

16 

7.20666 

0.0309309696 

0..400 

6 

20 

20.18198 

0.5534007984 

0.400 

8 

23 

23,12694 

0.5466209213 

0.400 

10 

1 

25 

22.04679 

0.3669595654 

0,500 

2 

8 

4.96920 

0.2391366628 

0.500 

4 

12 

13.04976 

0.6345655810 

0.500 

6 

15 

16.03409 

0.6202105433 

0.500 

8 

17 

25.16863 

0.9090016643 

0.500 

10 

19 

24.55579 

0.8243107298 

0.600 

2 

6 

13.11971 

0.9588261108 

0.600 

1 

-i 

9 

14.42509 

0. 892013025G 

0.600 

6 

11 

13.29194 

0. 7253265963 

0.600 

8 

13 

17.84469 

0.8364937782 

TABLE  5.2.3  (continued) 


p 

R 

D.  [■' 

0.600 

10 

14 

0.700 

2 

5 

0.700 

4 

7 

0.700 

6 

8 

0.700 

8 

10 

0.700 

10 

11 

0.800 

2 

7 

J 

0.800 

4 

5 

0.800 

6 

6 

0.800 

8 

7 

0.800 

10 

8 

0.900 

2 

2 

0.900 

4 

3 

0.900 

6 

3 

0.900 

8 

4 

0.900 

10 

4 

CM I -SQUARE 

SIGNIFICANCE 

13.95647 

0.5470410000 

6.42472 

0.7329438096 

11.10162 

0.8657526004 

9.98931 

0.7342226732 

4.82568 

0.0974859792 

14.96953 

0. 8161040822 

1.03482 

0.2071720727 

2.881S4 

0.2818034287 

1.48994 

0.0398398009 

4.S9512 

0.2907654804 

2.13612 

0.0234479888 

0.67942 

0. 2880246110 

0.23420 

0. 0281 120703 

4.99992 

0.8281970919 

2.59324 

0.3719796012 

2.34658 

0.3276973945 
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5,3  Random  Numbers  From  Continuous 
Distribution  Functions 

Sequences  of  random  numbers  from  the  exponential 

f(t)  =  c~X  ,  0  <  t  <  »  (5.3.10) 

and  normal  • 

f(x)  =  (1//27)  s'hx  ,  -«  <  x  <  »  (5.3. 2) 

densities  were  generated,  and  segmented  into  20  intervals  of  approximate¬ 
ly  equal  probability,  =  0.05.  Samples  of  m  =  1000  random  numbers  were 
generated  and  the  y?  statistic  computed  as 

20 

X2  =  £  (O.  -  m  •  p  )2/m  •  p  (5.3.3) 

i=l  1 

where  Cb  is  the  number  of  observed  occurrences  in  the  i*^  interval.  The 
significance  level  was  calculated  as  given  in  (5.2.3)  with  v  =  19.  The 
results  of  these  tests  are  given  in  tables  5.3.1  and  5.3.2. 

Sequences  of  random  numbers  from  the  gamma  density 

f(t)  =  (1/I'(n))  tn_1  e_t  ,  0  <  t  <  «,  (5.3.4) 

were  tested  in  a  somewhat  different  manner.  If  t  is  a  variable  with 
density  (5.5.4)  then 

y  =  2  •  -  (4  •  n  -  l)'"2  (5.3.5) 

should  have  an  approximate  standard  normal  density  when  n  is  large. 
Sequences  of  m  gamma  distributed  random  numbers  were  generated  and  using 
the  transformation  (5.3,5),  transformed  to  m  approximate])’  standard  normal 
variables.  As  before,  the  normal  distribution  was  partitioned  into  20 


RESULTS  OF  CHI-SQUARE  GOODNESS  OF  FIT  TESTS  ON  RANDOM 
NUMBERS  FROM  THE  EXPONENTIAL  DENSITY,  F (X)-EXP (-X) , 


D.F. 

CHI-SQUARE 

SIGNIFICANCE 

19 

20.76000 

0.6497876 

19 

11.24000 

0.0844989 

19 

17.64000 

0.4534074 

19 

22.68000 

0.7482694 

19 

13.68000 

0. 1979842 

19 

18.52000 

0.5120074 

19 

11.48000 

0.0933650 

19 

17.36000 

0.4345072 

19 

24.24000 

0.8128932 

19 

37.24000 

0.9925960 

TABLE  5.3.1 


i 
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RESULTS  OF  CHI-SQUARE  GOODNESS  OF  FIT  TESTS  ON  RANDOM 
NUMBERS  FROM  THE  NORMAL  DENSITY 


D.F. 

CHI-SQUARE 

SIGNIFICANCE 

19 

19.59699 

0.5808149783 

19 

25.42150 

0, 8528693863 

19 

18.29S36 

0.4974018037 

19 

6.71247 

0.0044039802 

19 

14.75770 

0.2621316817 

19 

19.91207 

0.6001134403 

19 

20.76446 

0.6500395713 

19 

23.50146 

0/7840244000 

19 

11.96570 

0.1129139927 

19 

18.0246S 

0.4792096579 

TABLE  5.3.2 


I 
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intervals  of  approximately  equal  probability  p^  =  0,05,  and  the  x2 
statistic  computed  as  before.  The  results  of  these  tests  are  recorded  in 
table  5. 3. 3. 


RESULTS  OF  CHI -SQUARE  GOODNESS  OF  FIT  TESTS  ON  TRANSFORMED 
GAMMA  RANDOM  NUMBERS 


SAMPLE  SIZE 

N 

D.F. 

CHI-SQUARE 

SIGNIFICANCE 

100 

50 

19 

12.301S6 

0.1276659399 

100 

100 

19 

14.27597 

0.2326276646 

100 

200 

19 

13.79290 

0. 20437271  OS 

200 

so 

19 

13.99520 

0.2160288470 

200 

100 

19 

20.25710 

0.6207407158 

200 

200 

19 

23.61244 

0. 78S558S803 

500 

50 

19 

23.13291 

0.7684615439 

500 

100 

19 

21.91446 

0.7114659110 

500 

200 

19 

32.02104 

0.9689153752 

TABLE 

5.3.3 

5.4  Discussion  of  Results. 

It  should  be  emphasized  that  no  test  or  series  of  tests  can  assure 

i 

the  suitability  of  a  given  sequence  of  random  numbers  for  a  particular 
problem.  When  possible  the  sequence  should  be  tested  on  a  similar 
problem  with  known  solution.  However,  the  results  presented  here  esta¬ 
blish  the  fundamental  reliability  of  the  generation  procedures  presented 
in  this  paper.  The  significance  levels  obtained  should  follow  the  uni¬ 
form  distribution  on  (0,1)  if  the  corresponding  null  hypothesis  tested 
are  valid.  The  above  results  arc  in  agreement  with  this  condition. 


CHAPTER  6 

DESCRIPTION  AND  USE 

OF  GENERATOR  PROGRAMS 

• 

6.1  GENO 

GENO  supplies  the  user  with  a  very  fast  procedure  for  the  generation 
of  independent  uniform  (0,1)  random  numbers  with  density 

f(u)  =  1,  0  <  u  <  1.  (6.1.1) 

The  method  used  was  suggested  by  Dr.  Rolf  Bargmann  [1]  and  is  a  form  of 
the  multiplicative  congruential  procedure, 

ur+1  =  a  •  un  mod  2 32  (6.1.2) 

where  a  =  513.  This  procedure  requires  that  Uq,  the  starting  value,  be 
an  odd  positive  integer. 

GENO  is  written  as  a  subroutine  with  entry  point  URAN.  It  is  called 
as  follows: 

U  =  URAN(IOD) ,  (6.1.3) 

where  "IOD"  is  the  starting  value  and  must  be  an  odd  positive  integer. 

The  desired  random  number  U  in  (6.1.3),  is  returned  to  the  calling  program 
in  single  precision  real  mode. 

GENO  requires  20  words  (80  bytes)  of  storage  space  and  an  average 
time  per  generation  of  25  uscc. 
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6.2  GEN1 

GEN1  provides  the  user  with  a  fast  procedure  for  the  generation  of 
exponentially  distributed  random  numbers 

f(x)  =eX,  0  <  x  <  “.  (6.2.1) 

The  method  used  is  that  described  in  4.1,  and  presented  by  Marsaglia  in 

[13]. 

GEN1  is  written  as  a  subroutine  with  entry  point  RANEXP.  It  is  called 
from  a  FORTRAN  main  program  as  follows: 

X  =  RANEXP (IOD) .  (6.2.2) 

The  parameter  "IOD"  primes  the  scheme  for  the  generation  of  uniform  (0,1) 
random  n ambers ,  "IOD"  must  be  an  odd  positive  integer  to  insure  the  proper 
generation  of  u.  The  required .random  number  x  in  (6.2.2)  is  returned  to 
uhe  calling  program  in  single  precision  real  mode. 

The  average  time  per  generation  is  approximately  70  usee.  Memory 
requirement  for  GEN1  is  307  words  (1228  bytes). 

6.3  GEN2 

GEN2  supplies  the  user  with  a  fast  procedure  for  the  generation  of 
ganuna  (a  =  n,  B  =  1) 

fi(y)  =  l/r(n)  •  ya_1  •  e“> ,  0  <  y  <  «  (6.3.1) 

distributed  random  numbers.  The  technique  used  is  described  in  4.2. 

GEN2  is  written  as  a  subroutine  with  entry  point  RANG AM,  and  is 
called  from  a  FORTRAN  main  program  as  follows: 
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Y  =  RAN GAM ( I OD , N )  .  (6.3.2) 

As  before,  "10D"  must  be  an  odd  positive  integer  and  n  =  "N”  in  (6.3.1). 

The  required  variate  is  returned  in  single  precision  real  mode. 

GEN2  requires  313  words  of  memory  space  and  approximately  65  •  N 
usee  per  generation  for  N  >  1. 

6.4  GEN3 

GEN3  supplies  the  user  with  a  fast  procedure  for  the  generation  of 
random  numbers  from  the  standard  normal  density 

f(x)  =  (l//2ir)  e~^X  ,  -»  <  X  <  ®.  (6.4.1) 

The  procedure  is  described  in  4.3  and  by  Marsaglia  in  [12]. 

GEN3  is  written,  as  a  subroutine  and  is  called  from  a  FORTRAN  main  pro¬ 
gram  as  follows: 

X  =  RANORM(IOD).  (6.4,2) 

The  parameter  "IOD"  must  be  an  odd  positive  integer,  and  the  required 
variate  X  in  (6.4.2)  is  returned  to  the  calling  program  in  single  precision 
real  mode. 

t  '  .  . 

■  Memory  requirement  for  GLiN'3  is  442  words.  The  average  time  per  genera¬ 
tion  is  approximately  65  uscc. 

I 

6.5  GEN4 

GEN4  provides  a  fast  procedure  for  the  generation  of  random  numbers 
from  the  binomial  distribution 

f(x)  =  (")  PX(1  -  P)n"\  0  <  p  <  1,  x=0.1f...>R.  (6.5.1) 
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GEKM  is  written  as  a  subroutine  with  two  entry  points,  BSHTUP  and  RANBI. 
The  user  calls  BSF.TUP  with  the  parameters  "P"  and  "N"  e.g. 

CALL  BSHTUP  (P,N) .  (6.5.2) 

BSHTUP  calculates  point  probabilities  using  the  recursion  relation 

f(0)  =  (1  -  P)N  .and  (6.4.3) 

f(x+l)  =  f(x)  •  P  •  (N-x)/((x+l)  •  (1-P)). 

The  four  sets  of  discrete  values  are  then  stored  using  the  procedure  de¬ 
scribed  in  Chapter  2.  BSETUP  need  be  called  only  once  for  a  given  "P" 
and  "N",  but  must  be  called  at  least  once  prior  to  calling  RANBI, 

RANBI  addresses  the  very  fast  scheme  that  generates  the  required 
variates.  It  is  called  as  follows: 

X  =  RANBI (10D)  (6.4.4) 

where  "IOD"  must  be  an  odd  positive  integer,  and  X  is  returned  in  single 
precision  real  mode. 

In  its  present  form  GEN4  requires  the  "P"  be  in  the  real  mode, 

,  0  <  P  <  1,  and  that  "N"  be  in  the  integer  mode  and  not  exceed  255.  The 
program  requires  948  words  of  memory  and  an  average  time  per  generation  of 
approximately  35-40  usee.  Figure  2.1  describes  the  flowchart  of  GEN4 
except  for  the  section  that  computes  point  probabilities. 

If  the  user  desires  random  numbers  returned  in  the  integer  mode,  he 
need  only  delete 


STMr 

199 

ST 

0, 

RLS 

200 

MV  I 

ri:s 

,  X 

201 

A  E 

0, 

RLS 

34 


and  change  the  designation  RANBI  to  IANB1  in  his  calling  statement  and 
in  statements 


170 

ENTRY 

RANK  I 

171 

USING 

RANBI 

,  15 

173 

DC 

CL7  * 

.  .RANBI 

174RANBI 

STM 

1  »  ^  , 

24(1*3)  • 

Storage  requirements  may  also  be  altered  simply  by  changing 
STMT 

146  C  7,  =  F ' 2000  * 

222  DS  2000XL1 . 

6.5  GEN5 

GEN5  supplies  a  fast  routine  for  the  generation  of  random  numbers  from 
the  poisson  density 

f(x)  =  AXe"X/x!  ,  0  <  A  <  »  ,  X  =  0,1, _  (6.5.1) 

Since  the  domain  of  x  is  infinite,  it  is  neccessary  to  truncate  the  distri¬ 
bution.  This  is  done  at  both  low  and  high  point  values  whenever  a  point 

-4 

probability  is  less  than  16 

GENS  is  written  as  a  subroutine  with  two  entry  points,  PSETUP  and 
RANPOI.  The  user  calls  PSETUP  with  the  parameter  A,  e.g. 

CALL  PSETUP (ALAM) .  (6.5.2) 

PSETUP  calculates  point  probabilities  using  the  recursion  relation 

f(0)  .  e'X 


f(x+l)  =  f(x)  •  A/ (x+ 1 ) , 


(6.5.3) 
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H 


and  stores  the  four  sets  of  discrete  values  by  means  of  the  technique 
described  in  Chapter  2".  PSETUP  need  be  called  only  once  for  a  given  X, 
but  must  be  called  at  least  once  before  the  initial  call  to  RANPOI. 

The  entry  point  RANPOI  addresses  the  very  fast  scheme  that  performs 
the  actual  generation  of  the  required  variate.  It  is  called  as  follows: 

X  =  RANPOI (IOD) ,  (6.5.4) 

where  "IOD"  must  be  an  odd  positive  integer  and  the  desired  random  number 
X  is  returned  in  single  precision  real  mode. 

GEN5  requires  972  words  of  memory  space  and  an  average  time  per  gener¬ 
ation  of  approximately  35-10  usee.  The  flow  of  GENS  follows  that  presented 
in  figure  2.1  except  for  the  section  that  computes  the  point  probabilities. 

If  the  user  desires  his  result  returned  in  the  integer  mode,  he 


need  only  delete 

STMT 

* 

190 

ST 

O.RES 

191 

MV  I 

RES, 

192 

AE 

0,RES 

and  change  the  designation  RANPOI  to  IANP0I  in  his  calling  statement  and 


in  statements  ' 

STMT 

161 

ENTRY 

RANPOI 

162 

USING 

RANrOI 

164 

DC 

CL?'  RANPOI’ 

165  RANPOI 

STM 

1,3,24(13) 

Storage  requirement  may  also  be  altered  by  changing 
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STMT 

137  C  7,  =  F'2000* 

211  A  DS  2000XL1 

6.6  GEN6 

GEN6  provides  a  fast  routine  for  generating  random  numbers  from  the 
negative  binomial  density 

Y  +  V-  ^  T  Y 

f(x)  =  C  x  )  P  (1-P)  ,  0  <  p  <  1,  x  =  0,1,  (6.6.1) 

and  r  >_  0. 

The  domain  of  x  is  infinite,  consequently,  the  distribution  is  truncated 
at  both  low  and  high  point  values  whenever  a  point  probability  is  less  than 
16-4. 

GEN6  is  written  as  a  subroutine  with  two  entry  points,  NBSETU  and 
RANE8I .  NBSETU  receives  the  parameters,  p  and  r,  and  calculates  point 
probabilities  using  the  recursion  relation 

f(0)  =  pr  (6.6.2) 

f(x+l)  =  f(x)  •  (1-p)  •  (x+r)/ (x+1) . 

♦ 

Having  computed  the  point  probabilities,  NBSETU  then  stores  the  four  sets 
of  discrete  values  using  the  procedure  described  in  Chapter  2.  NBSETU  is 
called  as  follows: 

CALL  NBSETU (P, I R) ,  (6.6.3) 

where  "IK"  must  be  in  the  integer  mode  and  0  <  "P"  <  1.  NBSETU  need  be 
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called  only  once  for  a  given  "P"  and  "IR",  but  must  be  called  at  least 
once  prior  to  the  initial  call  to  RANEBI. 

RANEBI  addresses  the  very  fast  scheme  for  the  generation  of  the 
required  variates,  and  is  called  as  follows: 

X  =  RANEBI  (I0D)  ,  (6.6.4) 

where  "IOD",  the  primer,  must  be  an  odd  positive  integer,  and  X  is  returned 
'in  single  precision  real.  mode. 

GEN6  requires  944  words  of  storage  space  and  an  average  time  per 
generation  in  the  35-40  usee  range.  The  flow-chart  of  GEN6  follows  that 
of  figure  2.1  except  for  the  section  of  GE.N’6  that  calculates  the  point 
probabilities. 


6.7  GEN7 

GEN7  provides  a  fast  procedure  for  generating  discrete  random  vari¬ 
ables  with  any  specified  probability  distribution  denoted  by  a  vector 
P  =  {Pj»P2>  . It  is  written  as  a  subroutine  with  two  entry  points, 
DSETUP  and  RAND IS.  DSETUP  is  called  with  3  parameters,  "P",  "IX”,  "N"; 
where  "P"  is  a  vector  of  probabilities  (real  mode),  "IX"  is  a  vector  of 
corresponding  discrete  values  (integer  mode),  and  "N"  is  the  number  of 
elements  in  vectors,  e.g. 

CALL  DSETUP (P,  I X,N’).  (6.7.1) 

The  discrete  values  of  the  "IX"  vector  are  stored  in  four  sets  based  on 
the  high  order  4  digits  of  the  elements  of  the  "P"  vector.  This  procedure 
is  given  in  Chapter  2  and  illustrated  in  figure  2.1.  DSETUP  need  he  called 
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only  once  for  a  particular  set  of  parameters,  however,  it  must  be  called 
at  least  once  prior  to  the  initial  call  to  RANDIS. 

RANDIS  addresses  the  very  fast  scheme  that  performs  the  actual  genera¬ 
tion  of  the  required  variates.  It  is  called  as  follows: 

X  =  RANDIS (10D) ,  (6.7.2) 

where  "IOD"  must  be  an  odd  positive  integer  and  "X",  the  generated  variate 
is  returned  in  single  precision  real  mode. 

In  its  present  form  GEN7  requires  that  the  elements  of  the  "IX"  vector 
be  in  the  integer  mode  and  that  no  value  exceed  255.  Also,  the  parameter 
"N",  the  number  of  elements  in  the  vectors,  may  not  exceed  256.  If  the 
user  wishes  his  discrete  variables,  members  of  the  "IX"  vector,  to  be  in 
the  real  mode  or  to  be  in  the  integer  mode  but  with  values  greater  than 
255,  lie  need  only  make  the  changes  indicated  by  table  6.1.  If  he  desires 
to  have  his  result  returned  in  the  integer  mode  he  need  only  delete  the 


following 

STMT 

- 

149 

ST 

0,RES 

150 

MV  I 

RES,  X 

151 

AC 

P,RES, 

and  change  the  designation  RANDIS  to  IANDIS  in  his  calling  statement  and 
in  the  following 
STMT 


120 

ENTRY 

RANDIS 

121 

USING 

RANDIS ,  15 

125 

DC 

CL 7’  RANDIS 

124  RANDIS 

STM 

1,3,24(13). 

£9 


Storage  requirement  may  also  be  altered  simply  by  changing 


STMT 

97 

C 

II 

168  A 

DS 

2000XL1 

In  its  present  form  G nN7  requires  904  words  of  memory  space  and 
35-40  usee  per  generation.  • 
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X 

X 

< 

* 

* 

< 

25 

o 

V 

*- 

■» 

a 

n 

2 

u. 

• 

• 

< 

— ■ 

T‘ 

U j 

cf 

»~ 

ir 

r-> 

* 

ft 

o 

• 

r» 

Uv 

i 

cc 

4 

< 

«— 

•m 

4 

3 

l^c 

C 

ft 

ft 

« 

*-» 

*_ 

u> 

# 

« 

tr. 

* 

«T 

torf 

to 

ST 

_ 

c 

ct 

C 

fc 

* 

ft 

«• 

-r 

X 

CL 

> 

f— 

c. 

< 

<r 

tr 

♦ 

» 

C 

rv* 

<r 

to 

X  00 

ISJ  ^ 

u 

-J 

• 

> 

-0 

i; 

H 

* 

z  z 

m 

• 

• 

ft 

D  < 

UJ 

» 

ft 

to  ft 

•l 

2 

c. 

c 

«<■ 

k- 

• 

-1 

cf 

r— < 

r» 

c  x 

«C  — J 

cn 

or 

V' 

af 

rr 

ft 

c: 

ft 

C 

u. 

« 

ft 

ft 

-J 

-J 

* 

»  UJ 

»  * 

« 

ft 

to 

• 

» 

af 

O 

r 

< 

UJ 

iz. 

■r 

LLJ 

« 

“ 

3 

X 

ft> 

r-< 

O  ft 

fNj  fO 

fF' 

o 

ft 

r 

ft  ft 

ft  CC 

* 

»/> 

3? 

cr 

CL 

c 

•Y 

UJ 

ft 

3 

•»* 

cf. 

UJ 

OJ 

< 

>  C3 

• 

C 

U. 

ti 

3 

3 

X 

X 

• 

-> 

».. 

4-> 

X 

o 

H 

t— 

1 

* 

»- 

— 

T 

a 

ft 

*— 

ft 

o 

I 

1 

1 

X 

l 

» 

Z  ^ 

ft 

u>  4- 

ft 

ft 

a 

ft 

> 

ft 

X  o 

> 

i 

1 

» 

•« 

1 

# 

u. 

3 

ft  c 

00 

to  ec 

X  V* 

to 

to 

£ 

< 

ft  cr 

ft 

X 

CO  o 

z 

X  — 

O  ■  Lt 

y  (t  a 

<  uj 

a  7  ft 

CC  -  c 

a*  y. 


ft  t".  •  2 


I 
.  u 
c  a  * 
<  <  < 
C  C  a 
ft  ft  tf> 


«t 

~  ^  X 
*  *  u 
~  ~  2 
wC  «* 


v/1 
!  O 


I  UJ 
UJ  S' 
Ux  ft 

a;  & 
D  a c 
C  => 

on  ft 

P  4 


-I  c 

a  c. 
x  r 
<  »- 
X  U- 

UJ  X 


c  ^  n  >t  «* 


3  a 
n  m  c 


*  <\i  m  ir>  «o  ft  a*  ?  C  ■—  <v. 

vrjru  fufvirgrvrurvru  m  <c,  , 


C 

O 


a 

o 

o 


CUC®!' 

^  -f  IT  C  ^ 

o  c  c  a  c 

O  C  c  o  c 

c-  c  o  o  c 


CD  CD 
•-* 

C.  3 


O  <>  C  O  rsi 

COi" 

a  c*  o  a  o 
cocoa 
C  C  c  o  c 


l/' 

o 


CJ  « 
*#■  • 


C  l 
%f 


o 

a 

c- 

c 


or  **  tc  u-  oOocd 

rw*J  *4  C-  *  -*  -4 

C  C  O  O  C  C/  c  c. 

U.U.UCUUUO 


r-  c. 

c  * 


a>  -c  o  o  < 

->*  o  x  c.  o 
ct  C  C  o 
a  u.  u. 

ft  r*  O  c  /“*  r-  r  n  ^  ,#i  n-  ^  •»  f.  r>  o 

—  c;  u  r\  r<  j  c  u.  r  o  f*1  (Ml 

C  <r  ft  ft  cr  c  r,  <  cc  ft  cr*  .x  a  uj  ft 

C  fA  O  U“  ^  C  If  c  f'  0  C  IT  C  IT  IT  «l 


I  C.'  O’  u  UJ  (\  -C  <3  U  f\  ~c  «IUjO>f 

ooc  ft  ft  C  --4  — *  — *  —4  fM  fX  r\,  rv  rn  ft  ft  ft  %*  s| 

w  t  C  t  u  C.  w  C  ^Ci  —  Cw  — 

C  C  C  ft  u  C.  "  Q  r-  O  C  1  C  O  C'  C  C-  l  O  T 

re  r  or^.  occe  our'ooc^C'r  o  ^ 

f  i  l  t  c  »  c  c  c  r  r  c  :i  r  o  c  c  n  ^  r 
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C 
u. 

c 

ft  u. 

9  Q 

or  w. _ 

ft  ow  »- 

<  o  tt  7 
a  w- 

•  a  r 

•  *  c 
*-  *-  r  c 

<  ^*zr 
►-  <■ 

v»  u  x  z 
o  •>* 
w  ft  C 
C_  Z  w 

•  ►*  or 

ft-  —  ft- 

A.  X  Q  V* 
ft.  u  3 
o  a  ft  r 
*u  3 

•  A  «  • 
«<  3  —  — 
O 


O  kw 


—  t  • 

Z  A  • 

T>  u-  • 

J  X  • 

r  o  • 

Co  • 

UO  # 

•  O  •  * 

<  * 

X.  —  • 
7<  J* 

—  #  Cl  • 

O  4  I  iT 

c.  c  v*  •  -o 
iw«  a  t  s 
a:  a  <  •  c 

—  X  •  -c 
or  w.  •  r 
kj  X  v  « 
vi»>  o  «  a 


£  o 


x  e- 
•  w  ft 


a  .a* 

«*  UJ  • 
'V*  <*  V*  • 

<  a*  o  • 


—  >  C  •  5 
-»4  -  ft  3 
e  x  4-  «  *- 
4  •  4  •  ^ 

»/•  4 

or  .Oft  • 

*  C  Z  •  u. 

xx  •  r 

ftwftC 
u  ft  X  •  C 
C  •  ft-  • 

—  ft  c 
ft-  C  X  •  at 
-ft  >  O  ft  Gk 
•  X  -  j# 

■  ft  •  -J  • 

X  C  ft 

« 

—  <  •  — 


u,  •  ft 


V  C  •  w 


O  ft 

o  « 

4  w 
or  A 

Z  AC. 
A 

?  «# 

• 

u.  — 
O  r 

VI 


VI 

O  ft 

*-  a? 

c.  — 


o 

X 

o 


X 

c 


D  . 
o  : 


_#  "  o 
«  •  X 

►-  r  ^  < 

C  *•  7  a 
W  »  X  u. 

icS: 

ft  u.  «  or  • 
ac  a  or  o  • 
C  *•"  <r  v  • 
«  fc 

ft-  X  T  ►*  c  o 

•  CCft  X 

—  c;  •  z  c 

or  X  A  —  A 

ft  A.  A»  ft  3 

U.  c  c  o  -  ■ 

—  x  —  —  a:  a 

O  r-  O  C  C  O 


AJ 

X  C 

O  — 

r  c 
•»  »* 
a  ft 
*  a 


—  > 

or 

•  x 

14.  tf  O 
A  Z  C 
ft  / 

®  X  ft 

-  c  a 
c  X 

Aft  ft  14.* 
n  a  ►- 
I  C  Uj 

—  or  • 
C  A  O  • 

.  —  ft  ft 
ft  —  ft- 
L  I  O  6 
V*  A  ft 

ft  o  ^  r 
K.  C  —  ft 
—  A  ft  C 

<  c  »-  « 

«r  —  «.  ft 
ft  O  O  or 


v.  ft  u. 
>  v*  —  .j 
uj  or  or 

or  ft  ft 

e  >  *•• 


ft  ^  UA 


a  o  ft* 
C  v»  u.- 
«  ~  at 

.  aCy 

>  z  z  - 
~  -  c 

.  <  < 

*  33  tf  ft 
O  C  or 


C 

ft  09  | 
k  ft 

—  aft  • 
o*u.  P* 

O  J1 
ft)  ft  u.  • 

~c  O  ♦ 

«*»  K 
UJ  7  ft 
C  3*^ 
a*  Z  Q  3 
C  ft  —  A 

—  ft  u; 
*>s  «  a 
X  u. 

—  >  ft 

♦  a  z  z 
►  ft  C  ~ 

.  UJ  —  vj  ft 

4- 

cd  a  7  t 
ft  «0  O  C  a 

-  «  w  c  « 

or  *■ 

ft  or  ft  c  * 

>«L 

o  C  a 
u;  or  u 

*-  O  C 
tw  u.  . 
or  C  C  »- 
ujj2 
ft  uj  — 

—  X  *S.  C 
C  c*  a 

Z  ft 

«  •  »*o 

ft  ft*  —  * 
or  at  ®  — 


O 

u;  Z  C 
a  —  ft 
•  ft 
lul  3 
or 

O  O  u* 
h  t-  or 
ft  O 
lM  X  ft 
a  u  ft 
2  u. 
ft  Ot 


—  O 
o  c 
ft 
c 

X  UJ 

«  — 
ft  w 

v  o 


ft  z 
X 

S* 

ft  ft 
C  a 
X  ft 

z  — 


Sd 


z  c  ■ 
c  > 


-J  « 

■  r*  -v  ^  * 

«CC  «  i 
.  X  3  « 


4  &  4. 

VJ  ~  V 

•  ft  u. 
W  7 

•  C  C. 

— *  *-  c 


I  I 
I.  I 


.'-7 

ft  —  ft 

ci  ft  a  i 

*  ft  4  : 
*•  U.  X  I 
■  I 

I  I 

•  u. 

I  -J  4 
UJ  &  ( 

o  x  : 


/  2 

•  ft  <  _ 
ft  c*  * 

>•  t 


Z  A  t-  V 
u.  D  O  ( 


Z—  »w  -C  —  — •  — 

-ft—  A  7  ft  ft  ft  A  ft  — 

r*  a  4  *-  ft  <7  4  4  ft  ft  —  •  4  —  — 

c  N  ft  .ftrs  k  -c—  u  -u-v  V  V  vw\  4  ft  o  ft 

-••A.  n.  ft  uZ*  w  UZ»  uj  ll  7  <4  —  o  uj  -  uj  *  —  ft  —  Z  — 

—  —  *Cft  -  ft  ■  •  C  ft  •  ft  •  »  lOC^cnraD  *  .u  ft  ft  a  -  r*  j  a-  <c 

-J  ~J  •  U.  »*.-.*  »  W*  *  •  —  -  -  A  *•—  -  -  •  «*u,»  -V.  -  •••—  4  •  • 

*  U*  —  ►>  C  ft  A  ft  ft  ft  ft  — ftA*uft  ft— ftft  —  ftft—  ft  ftftftir.cro  —  —  AAjftft  —  ftftA 


a-  X 

K  U 

U.  z 

Z  ft 
ft  or 
or  « 


•ftAftr-ccC— ruftftftft^-i 


S 

ft 

— 

« 

8 


r  o 
o  o 
c  c 


J  o  -I  O  C  —  Of 

J  U  u  UU  JU  u  O  o  *  or*  >ll<Z  VJ  4  B-  JO  K  (t  J 

'  OV'OC-^KiAOC^-J/ur'tft  J  £  JWVMJCM^4U 


«  «  »«IL  z  z 

>*-r  90— AftftAftAftcC  — A.ftftftftr-«9C  —  Kft.  ft  v 
.ruAJAftftftftftftr  ft  J-  JA  ft  ft  ft  ft  ftftft  ft  ftftuMA«rv>AA 


c 

o 


c 

c 


ft  N 
U.  U 


C  ft  1  ft 
C  u  V  U. 


e  o-  r  a  r  a 
C  —  u.  o  ft  A 


C  ft  O  ft  ft  A  A  —  • 


A  U. 
ft  C* 

ul  c  u  u 


«  C 

t  C 


ft  c 

c  o 


-»  uo 
O  o  o  — 
C  .0  3 
C  ft  5  ft 
C  C  C  O 
t  C  C  r* 


NC  4  O 


c  c 

w  C 


AftAuftoftftftOlUAO  O  ft  ft  ft  «  «ft 

C.ACO  -|CCI'CUNC<  AC 

-  C  UC.AOC  JftCftft  c  —  —  —  G  —  O 
CUftUl4U.ftWCWU.UU.  ft  W  U  W  C*  U  c  H  .  » 

'3'“‘'^*ftftft',,,COftftftCft—J-ft^ftCftftftAA*** 
Al»ftftft4>AWft  ftft  ftAftAAftO**aU»Ars,A'rceftftft 
9ft  Xft«.u>C^4ft«C  f-ft  —  orAfta.ftA 

«  ftCAV^rtfftC-A^A^—  «v9ft<7CAft— Aft—  —  A 


•  •  X 

C  Ml  O  • 

*-  w  T  O.  « 

k-  0  «  J  < 


w  V*  cf 

z  u.  *-  > 

M.  •-  >  <  ft  i 

o  >  » 

cm.  •  X  : 

►-  *-  o  o  • 

•  K  Cui  7  < 

Du.  htf  JV<( 

Z  uw  2*a 

X  XV*  —  tf*  V 

w  «r  •  —  —  c  < 

z  c  e  »-  vi  u. 

U-  O'  k~  U.  c  V 

Z  ^  Z  X  v  -*  * 

k»  W  If  V  •-  o  <  -#«. 

k>Ca<7k>4au.C 
Ck-^CCatk«C 
CV>UUC«8  < 


Ik.  r  u.  <  P"  -J 

»-za  <  -  u.  ■ 

•  U  »c  JlNC  •  . 


X  -V*  X  »-  -*  O 

SJ  •  •  04  ft  I  4 

i  /  -•  4  ra  •  a 

I  *1  "»  O  «i  C  P*  NO  *K  ( 

■  a  •  •  •  d  »  •  *  •  * 


M  UJU_*-ctaou  >  *-  «  tc^u  Ml 

mi  cuaivi^iAHCjij^virxjiuczvi, 


o  *•  z 


ftyffto  «M**o«cceOtt  u  <«« 

■  r  t  ^  oN^i'-cr.  m 

cno--  r.  -r.  oc**r*-  —  c  — 

eoccc  coooooc  t*  o  o  o  c  c 

coccc  oc^scococc  eeo 


M- 

4  ♦  c  c  c  Nuu*/cocro  *  m.  o 

ft  ^  C  M  N  C<t}tf*KCOCNC^  |ku»k 

C  r»  C  *«  -«  orc-ccj-c*-  mC- 

ft  ft  ft  ft  ft  Oftftftft-*7v  ftft  ft  ft  ft 

o  ^  r  c  r  4rk*r*  ?k*f*r.c  ccn^co^mt' 

n  «ir>  c  4  ft  « 

►  /MjOfc«|k^C\r«Ck-0  i  «*v  »>ur  c 


ft«Ut»4ftftft-N44C4ftUO-l  Ut  4  < 

M  C  »  ft  ft  ft  ft  <  T  C  r  CWU  WC  C  o  c  C  u.  w  W 

o  c  o  c  t  t  t  «-;•  2  ^  c.  o «-  c  c-  i  c  c  c  c  c 

coccc  o  v  c  c  o  p  c  n  o  or.  c  00  1  o o 

tooocccccccCcr  r  o  o  30  i c  c  o 

ococcooccnr.  rnc  r  r.  ^  »  c.  c  rrc 
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tcc  object  cooe  *00*1  ion**  stmt  source  statement 
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» 


*  —  — 

ir— —  c 

-*fN*«vcr*(r«A,c 

cv  f*“  ev  rev  «  a  r- 

c«KON<?ET 

»r  x  a-  ■*>  ®  c  c  e 
co-eco-eeo* 


®  c 
c  C 
r  «r 
a.  r: 

v'  c 

c  — 
•#  *v 
g  g 


of>c  r  c  do  c  cc 


c*  p  ^  ^ 

r~rvE»‘*/^R»f'^C— •OC7r<.^<*>w'tf  V 
*/  4*  O  A-  <A 

r»f>*£r..A  —  a:~«.c  *v  ®  *  a.  •*  xycvi#iW'«AK.^*« 
*#  A  a.  f  N  /  N  C  "  Chew'S.*". 
/'Tr^C'OR.  <eJomririo-<rN^*fC«‘r'0 
w*  e 

iNC^^OfT.  «f  C  I/'  C  >/•  C-  ^  ? 

*^Arr?«\r^irfKNj.  Po  C  ^  en  r»  »e  r 

s*  *r  g  >s  s>  <  -c  **  ■*  -e  4*  c  c  v  <oaa.*»a.a.a  ►*  a- 

c  c*  e  e  c  a  o  c  e-  c*  c  e  pff  if  ^ 

rrr»^ooccoooopoo#  c  d  f  c  r  <*•  c  o 


TE.tftO 

—  a  c* 
r.  — 

IT-  —  tv 
<  •# 
N  1E  /> 

4  *  \ 

«j  *r  ir* 

N  K  K 

G  G  G 
•  »  « 
r*  r»  r* 


o  tv  o  * 

»•  y  k  «  n 

«•  C  4*  ~ 

cc.  4  ir  it*  — 

•w*  *•  •-•  c  *r 

r  -  ^ 
m  a  n  ir 

5  ^  K  C  <•* 

»  O  <  N  E* 

E  N  S  ^  K 

C  C  O  G  G 

•  *•11 
OC  CTO 


Run  *M 
MCTm. 


^  «r  il' 
r  *-  » 

®  e  — 


G  G  G 
•  *,* 
cc  b 


c  c  n* 
*  c  •* 

«/>  40  < 

a.  e  a- 

■T  N  < 

•C  G  G 

r-  *  g 
■r  ^  g 
«■  r  or 

A  A  ►- 


c  c  a 
•  •  • 
000 


<rc 

•#  c  c  dC 

TA  O.  40  TA  A*  ® 
\Q««N/ 
C7-  if  OO-eN 
C  N  <  K  O 
f~  h*  ~  •*  G  S3 
TV  ^  K  (f  »-  C 
^  C  C  (*  CC 
T»  f.  K  T*  a  * 
G  G  G  G  G 
•  •»••• 
c  c  c  o  r  o 


eocc*oc*cccoec.ccoococcooecco  0000c  c  c  <1  CA'Cc.ceoecccc.oceoncocc 


WIT  IW  urIL  U  Ul  W  II  ILUJIUWU.btUU.UU.  U.  Uj  14/  U<  U.  U'UMuIN  WU  U.  U-  u>  U  U  CUMi/U'IU  kUlLU  U-'IL  II  U.U.  IlIU  (U  U'U'  tiUf 
UMf  »E  IT  VIE  Irt  4Mf  N"  ✓  VI  ^  V4  VU"  4f»  «A  V"  TMfl^lTlf  Nl  *f  lE  /  If  4/!lC*04rt»040*f»%Air4f*  If  4T>  lA  VHf  ^  (/*«/)  if-  tflVlV>V4tAt/l 

*4t**<*<*«**-i**t^**-*<«tf<<<*<«’<*««r<<<«*<<«**<*<<<<«<<*f**<*-T«*< 


«C  4J  -CvJ-Ouww.^  V  -  WNW<  ^4Ku.0a.ff«.M«4K0i*iiA<C0*E-Ctt<UWNlAU<E* 

NftCT*-rOif<«<^:E*v0l<«wf'\Kuj^E-<f>if^ceci-^».ifirN4('CUlM3if«CU.4<*'ri0'O 

o  m  -  N  if  e.  *r  x  i  if  /  ?  c.  .*  cd  *»  c.it  r  -  c  f  *  *.'  :  r>  j  -f  f  f  if  e  v  t  j  NNM'.i.-'1i^1?  roM 

«|kC«w^*M.*C  -r  f/M  "f  r  r«  O  k.  —  X  w  C.'  C  c-.  3  t:  C‘  —  TV  <C  CC  -*  <  E*  <  IE  t;  iu  «  C  r*  •)  U  o  4  C  4  j'i  ^  o 

M  U-CE  E>  <  *;r  *  i‘  T  '•*  %  v  l»  QC  ™  »  T  I.-I  «  C  T.  E*>  f  <i:<C4’'/J<TCT<CE'C<ttOUUUU.wOtfffl 

W.  »A  «C  ^  U#  C-  *-  »\.  K  X  E  A  V  f  ff.  Jlf<iliN<(M,W**'TE.ffeC  ~  S'  c  u  U  C.  V-  r  EM  T~  «f  lO  <  A  «<f<XOOU-U. 
—  A#  A.  TV*  IN  *r.  fiRif  |M  ^  «r  <  If  S  S'  <3^5  «NNKNK  «T  *r  X  f  C  r  s  C  c  ff  o  c  C*  ?  c  «  4  •(  <  4  4  < 

ILklLku-wUlLlLUkUlL^bU  UWVUU.UUlLU.ILkU  U  11  ^  ^  w  k  U.UU.«i.U.U.U.U.UU.U.U.U,U.U.U.U.U.U.U.U 
XKMWMMKMKKHKWMMKMITKMMKKMiTKWWKMWMKMKKKMKMKkTKKKKXKMKMKKXW 
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■  -f  « 
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wi  tr.  «> 

• 

a  ■—  r-  C 

wax 

Z  X 

Uj  —  —  -1 

►*  < 

•-  c 

<  tL  C  +- 

?  Im  c 

w- 

•  o-  —  X 

VIUC 

1  4 

u#  •  •  * n 

X  X  U. 
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